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ABSTRACT 


The high speed associated with flight at very high altitudes 
introduces many new problems to the missile designer. Not the 
least of these is overheating of the missile structure due to fric- 
tional heating. 

At the extreme altitudes encountered during flight of sounding 
rockets or missiles, the atmosphere can no longer be considered 
as a continua, and account must be taken of the individual mo- 
tions of the molecules comprising the atmosphere. 

Using the methods of kinetic theory, calculations have been 
made of the temperature of uncooled flat plates traveling at high 
speed in the upper atmosphere. The calculation may be ex- 
tended to bodies of arbitrary shape by considering them to be 
comprised of a number of flat plates. 

It is pointed out that an effective method of cooling is to ensure 
thermal contact between portions of the body inclined at positive 
and negative angles with respect to the flight path. 

The effect of solar radiation on body temperatures is shown to 
be increasingly important as the flight altitude increases. 


INTRODUCTION 


Mss OR AIRCRAFT FLIGHT at extremely high 
altitudes introduces many new problems to the 
designer. One of the more important problems that 
arises is caused by the skin temperatures that are at- 
tained at very high velocities. The severity of the prob- 
lem is commonly exemplified by consideration of the 
disintegration by combustion of most of the meteorites 
that enter the earth’s atmosphere at high speed. 

In reference 1, Tsien has summarized the meager 
amount of research that has been done in the field of 
tarified gas kinetics as applied to aeronautics. It is 
shown that the general field of aeronautics may be 
divided into three major régimes and that the division 
lines, although poorly defined at present, are marked 


Presented at the Joint Session with American Physical Society, 
Sixteenth Annual Meeting, I.A.S., New York, January 29, 1948. 
* Aeronautical Research Scientist. 
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Heat Transfer to Bodies Traveling at High 
Speed in the Upper Atmosphere 


JACKSON R. STALDER* anp DAVID JUKOFF? 
Ames Aeronautical Laboratory, N.A.C.A. 


by the ratio of the mean free molecular path length to a 
characteristic body dimension. 

The three main fields of aeronautics, or gas kinetics, 
may be described briefly as: (1) the region of conven- 
tional aerodynamics in which the mean free path length 
is negligible in comparison with the body size; (2) the 
region of free-molecule phenomena in which the mean 
free path is large with respect to the body size; and 
(3) the so-called ‘“‘slip-flow’’ régime, intermediate be- 
tween the two prior régimes. 

In conventional aerodynamic theory, the first basic 
assumption that is made is that the fluid under con- 
sideration may be treated as a homogeneous, continu- 
ous medium. This assumption certainly should be- 
come invalid at very high altitudes because of the ex- 
treme rarefaction of the atmosphere. A measure of the 
atmospheric density is furnished by the average length 
of path (mean free molecular path) that is traveled by 
a molecule between successive collisions. The mean 
free molecular path has been computed by Johnson 
and Slack? for altitudes up to 150 miles. The result of 
their computations is shown in Fig. 1. The existence of 
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a large mean free path should not be construed as an 
‘indication that the number of molecules per unit vol- 
ume is small, in the ordinary sense of the word. For 
example, at a reduced density at which the mean free 
path is 10 ft., there are still 10'* molecules contained in a 
cubic inch. 


The work of this report will be confined to a study of 
the heat-transfer processes occurring in the region of 
free molecule flow. The lower limit of the ratio of mean 
free path to a characteristic body dimension at which 
the derived theory becomes invalid cannot be stated 
at this time for lack of experimental evidence. The 
numerical computations, however, will be confined to 
altitudes of 75 miles and higher. The estimated mean 
free path at 75 miles altitude is 1 ft. 


The specific purpose of the present report was to 
develop a method whereby the surface temperatures 
of uncooled bodies traveling at any speed in a rarified 
gas could be determined. Although this objective was 
realized, it should be pointed out that, to the authors’ 
knowledge, no experimental work has been done in 
high-speed rarified gas streams. Consequently, the 
numerical values obtained should be regarded with 
some reservation until experimental evidence is ac- 
cumulated concerning the validity of the assumptions 
that are made in the analysis. 


NOMENCLATURE 


A = constant of integration in Maxwell’s equation 

E = total absolute velocity, ft. per sec. 

Cx = total absolute velocity in X direction, ft. per sec. 

Cy = total absolute velocity in Y direction, ft. per sec. 

Cz = total absolute velocity in Z direction, ft. per sec. 

e = base of natural logarithms, 2.7183 

E; = energy of incident molecules, ft.lbs./ft.?-sec. 

E, = energy of molecules at plate temperature, ft.lbs. /ft.?- 
sec. 

E, = energy of re-emitted molecules, ft.lbs. sq./ft.?-sec. 

Er = energy of rotation, ft.lbs./molecule 

E, = translational energy of molecules striking front of 
plate, ft.lbs./ft.?-sec. 

FE; = translational energy of molecules striking rear of plate, 
ft.lbs. /ft.?-sec. 

E, = vibrational energy, ft.lbs./molecule 

f(V) = Maxwell’s distribution function 

g = gravitational acceleration, assumed constant at 32.2 
ft./sec.? 

4) = number of degrees of freedom of motion 

Jo = normal solar radiation at top of atmosphere, 93.4 
ft.lbs./ft.?-sec. 

k = Boltzmann constant, 5.66 X 10°™ ft.lbs./°F. per 
molecule 

m = mass of one molecule, slugs/molecule 

M = molecular weight, lbs. /lb.-mole 

My, = flight Mach Number 

n = number of molecules striking front of plate per sq.ft. 
per sec. 

n’' = number of molecules striking rear of plate per sq.ft. 
per sec. 

N = number of molecules per unit volume of gas 


No = Avogadro’s number, 2.73 X 10% molecules/lb.-mole 
Heat removed from body by an internal cooling sys- 
tem, ft.lbs./ft.?-sec. 


S 
ll 


R = universal gas constant, 1,544 ft.lbs./lb.-mole, °p 
absolute 
R; = total incident radiant energy, ft.lbs./ft.?-sec. 
R, total re-emitted radiant energy, ft.lbs./ft.?-sec. 
43 = surface area, sq.ft. 
i = absolute temperature, °F. absolute 
T; = temperature of incident molecules, °F. absolute 
T, = temperature of body skin, °F. absolute 
T, = temperature of re-emitted molecules, °F. absolute 
LU’ = absolute mass velocity, ft. per sec. 
l’y = component of absolute mass velocity in X direction, ft 
per sec. 
ly = component of absolute mass velocity in V direction, 
ft. per sec. 
Uz = component of absolute mass velocity in Z direction, ft, 
per sec. 
V = absolute velocity due to thermal agitation, ft. per 
sec. 
V,, = most probable molecular, speed, ft. per sec. 
Vy = component of thermal velocity in X direction 
V,; = component of thermal velocity in Y direction 
Vz = component of thermal velocity in Z direction 
space coordinate 
y = space coordinate 
A = space coordinate 
a = accommodation coefficient, dimensionless 
B = reciprocal of most probable molecular speed, sec. per 
ft. 
Y = ratio of specific heats, dimensionless 
€ = emissivity, dimensionless 
6 = plate angle, deg. 
v = number of adsorbed molecules per sq.ft. per sec. 
o = Stefan-Boltzmann constant, 3.74 107" ft.lbs. /ft.*- 
sec.-°F.* absolute 
T = time, sec. 
¢(a) = probability integral, (2/\/7) dX 
x = dimensionless group defined by Eq. (7) 
x’ = dimensionless group defined by Eq. (9) 
y = dimensionless group defined by Eq. (14) 
vy’ = dimensionless group defined by Eq. (19) 


ANALYSIS 


In order to determine the skin temperature of a body 
moving at constant speed through a rarified gas, it is 
necessary, fundamentally, to make an energy balance 
on a section of the skin. Energy may be added to, or 
subtracted from, the body by three distinct processes: 
(a) molecular energy transport to and from the body; 
(b) radiant energy transport to and from the body; (¢) 
energy added to or removed from the body by proc- 
esses occurring within the body. Thus, we may 
write a simple energy balance on a section of the body 


as: 
E.+R,=£,+R+0 (1) 

Each term in this equation may be considered separ- 

ately. 

Molecular Energy Transport to the Body 


It is a matter of common knowledge that any body 
surrounded by a gas is constantly bombarded by mol- 
ecules of the gas. It is this transport of momentum. 
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infact, that causes a pressure to be exerted on a body or 
on the walls of a container in which the gas is con- 
fned. Likewise, the property that we call ‘‘tempera- 
ture’’ is merely a manifestation of the magnitude of 
the mean velocity of the molecules. It is evident that 
the molecules possess, by virtue of their velocity, a cer- 
tain kinetic energy. Thus, a body immersed in a gas is 
constantly receiving and emitting energy from the 
bombardment of impinging molecules and the re-emis- 
sion of the same molecules from the body. If it were 
considered that the molecules carried no energy ex- 
cept that due to translational velocity, then, we should, 
by examination of their incident and re-emitted veloc- 
ities, be able to determine the net increase or decrease 
of the kinetic energy of the re-emitted gas and, hence, 
determine the molecular energy added to or sub- 
tracted from the body. In the ordinary case, however, 
the impacts of the molecules become obscured by the 
collisions of the molecules among themselves, and the 
problem becomes exceedingly complex if treated on a 
microscopic basis. On the other hand, if we consider 
a body immersed in a gas that is of low enough density 
so that collisions among molecules in the vicinity of the 
body are rare in comparison with collisions with the 
body, the problem becomes susceptible to analysis by 
the conventional methods of kinetic theory. 

If it is assumed that the gas through which the body 
is moving is so rarified that the motion of the impinging 
molecules is essentially unaltered by collisions with re- 
flected molecules, then it follows that the velocity dis- 


gas in thermal equilibrium—the Maxwell distribution. 
The Maxwell distribution law for a gas at rest states 
that, of a total of N molecules per unit volume, the 
fraction that has velocity components lying between the 
values Vy to Vy + dVx, Vy to Vy + dVy, Vz to Vz + 
Vz is given by 


{(VxVyVz)dVxdVydVz = Ae DX 

dVxdVydVz (2) 
where Vx, Vy, Vz are the velocities of thermal agita- 
tion with respect to a coordinate system, X, Y, Z, fixed 
in space. 

To determine the energy transported to the body by 

the impinging molecules, we must first determine the 
iumber of molecules that strike the body. This is done 
as follows: 
Assume that the body is stationary and that the gas 
has a macroscopic velocity, U, with components Ux, 
ly, Uz with respect to the fixed coordinates. Assume 
the body to be a flat plate, and fix the plate at the 
origin of X, Y, Z so that the X-axis is normal to 
the plate. The mass velocity U makes an angle @ with 
the plate. Fig. 2 is a diagram of the coordinate sys- 
tem. 

The absolute velocities may be written as the sum of 
the thermal and mass velocities; thus, 
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tribution of the impinging molecules will be that of a. 
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x 
Fic. 2. Coordinate system used in analysis. 
Cy = Uxt+ Vx, GQ=Uyt+ Vy, GQ = Uz 
Consider a small section dS of the plate. All of the 


molecules that strike dS in 1 sec., that have velocities 
lying in the range Cy to Cy + dCy, Cy to Cy + dCy, Cz 
to C, + dCyz, will lie in a cylinder having dS as a base 
and with a length in the direction of Cx, Cy, Cz of C, 
where 


C = + Cy? + Cz! 


The volume of the cylinder per unit time is CydS. The 
number of molecules of this kind, contained in the 
cylinder is then 


where f’(C) is Maxwell’s distribution function for the 
total velocities. Then, to find the total number of 
molecules striking dS per sec., we must integrate 


ff’ (C)NCydSdCydCydCz, 


between the limits Cy = 0 to ~, Cy = —~ to , Cz = 
—« to ». Molecules having a negative absolute 
velocity in the X direction cannot strike the front side 
of the plate. Then, 


n=NS2.S. Se” f'(C) CxdS dCxdCydCz (3) 


The distribution function for the thermal velocities 
may be written as 


f(V) = Aen Ux)*+ Uv) 
Now, by definition (cf. page 45, reference 3), f(V) = 
f(C) and dCy = dV x, etc. Thus for a unit area, 


CydCxdCydC, (4) 


This equation can be integrated, the result is 


—p?Ux? 
= 


2/78 ©) 


+Su+ 
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Fic. 3. The dimensionless group, x, as a function of (U/V»,) X 
sin 4 


The factor 8 may be shown to be related to the most 
probable molecular speed V, as 6 = 1 / Vm, While V» is 
given by V,, 2gRT. A is related to B by A = 
83/n'*. Also ‘Us may be written as U sin 6. Eq. (5) 
may then be written as 


sin? - U 
vs 


sin 0 +o sin (6) 


This is the result obtained by Sanger‘ and rederived 
by Tsien.' The factor U/V,, is related to the flight 
Mach Number Mp) by 


% = 


U/Vm = V(y/2)Mo 


A plot of the dimensionless group, x, defined as 
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Fic. 4. The dimensionless group, x’, as a function of (U’/Vm) X 
sin 0. 
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sin 6 [ +o = sin ») | (7) 


To determine the number of molecules striking the 
rear side of the plate, n, Eq. (4) is again integrated, 
except that the limits of Cy are changed to Cy = 0 to 
The result is 


is shown in Fig. 3. 


—o, 


sin sin | (8) 


The group, 


in? - U 
m 


is plotted in Fig. 4. 


To determine the kinetic energy incident upon the 
plate due to translation, Sanger, as shown by Tsien, 
multiplies the mass of gas striking the plate, nm, by 
the sum of the kinetic energy due to the mass motion 
of the gas, U?/2, and the average kinetic energy due to 
translation, which may be shown to have the value 
8/, RT. This procedure, however, is an oversimplifi- 
cation of the problem, since introduction of the term, 
%/, RT, presupposes that those molecules that strike 
the plate contain, on the average, the same energy as 
those in the gas considered as a whole. A more exact 
calculation of the incident translational energy may be 
made as follows: 


Referring again to Fig. 2, the kinetic energy per 
molecule contained in a aan placed on dS is mC?/2. 
The kinetic energy of those molecules in the cylinder 
with velocity components between Cy to Cy + dCy, 
Cy to Cy + dCy, Cz to Cz + dC; is 


(m/2)NC?Cxf(C)dSdC (10) 


The total kinetic energy of the molecules striking the’ 


front side of the plate is obtained by integrating Eq. 
(10) over the limits Cy = 0 to », Cy = 
Cz = — ~too,. Then, fora unit area, 


= (m/2)N S 2 Jo” C?Cxf(C)dC 


— o to ~, 


(11) 


However, 


C? = Cy? + Cy? + 


so Eq. (9) may be written as. 
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m Ty)? >—~[y)2 — Iz)? 
E, = Cx(Cx* + + Cz2)e7 Ux) +(Cy—Uy)*+(Cz—Uz) 1d CydCydCz (12) 
The result of a tedious integration is 


1 1+ 3/2 Va BUx[1 + 6(8Ux) | 
Be = 1 = v2 
nl tes + VrBUx[1 + 


+ (13) 


where » is given by Eq. (6). Using the relations, 8? = m/2kT, 8 = 1/V,,, Ux = U sin @and defining a dimension- 
less group as ; 


U 
Vr sin of + sin 6) sin @ 


L+ sin of + sin 9) 


m 


Eq. (13) may be written as 
E, = n[{(mU?/2) + pkT] (15) 


When the mass velocity U is taken as zero, Eq. (15) 
reduces to 
E, = 2nkT (16) 


which is the usual result obtained for the translatory 
kinetic energy crossing an imaginary plane drawn in a 
static gas. When U is large, the equation reduces to 


E, = n[(mU?/2) + 5/skT) (17) 


A similar computation yields, for the incident trans- 
lational energy on a unit area of the rear side of the 


plate, 
= n'[(mU?/2) + kT] (18) 
where 
U ( sin? @ 
sin 6 E sin 6) sin 
= 1 +—— 7 (19) 
1— Vez sin 6 sin 0) 
42.0 
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Fic. 5. The dimensionless group, ¥, as a function of (U/Vm) X Fic. 6. The dimensionless group, ¥’, as a function of (U/Vm) X 
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The value of E,’, of course, reduces to 
E,’ = 2n'kT = 2nkT (20) 


for U equal to zero and to E,’ = 0 for large values of 
U. The dimensionless groups y and y’ are shown in 
Figs. 5 and 6, plotted as functions of (U/V,,) sin 6. 

If the gas is considered to be composed of hard- 
sphere molecules, then, for a monatomic gas, the total 
kinetic energy incident upon the plate is given by Eqs. 
(15) and (18) for the front and rear sides, respectively. 
However, if the gas is composed of diatomic dumbbell 
molecules, each molecule may carry an additional 
amount of energy because of rotation and vibration. 
By the principle of the equipartition of energy, the 
maximum amount of energy carried by a molecule is 
(j/2)kT where j is the number of degrees of freedom of 
motion. A rigid dumbbell molecule may rotate about 
two mutually perpendicular axes; hence, the maxi- 
mum kinetic energy of rotation per molecule is 


E, = kT (21) 


It is conceivable that collisions between dumbbell 
molecules could cause the atoms composing the mole- 
cule to vibrate or oscillate along the line joining their 
centers. In this case, the total energy of vibration, 
composed partly of kinetic and partly of potential 
energy, would have the value 


Some question may arise as to the correct procedure 
to follow when considering a gas composed of two or 
more constituents. It may be shown that the results 
that have been obtained for a homogeneous gas may be 
extended to mixtures of gas by calculating m and EF, 
separately for each constituent. This is a result of the 
fact that in a mixture of gases in equilibrium, each gas 
has the same velocity distribution that it would have 
if the other constituents were not present. 

The result obtained in Eq. (15) shows that, in effect, 
the front side of the plate selectively collides with 
molecules possessing higher-than-average velocities. 
It is pertinent to inquire if this does not mean, also, 
that the faster moving molecules carry more-than- 
average rotational and vibrational energies. It is 
pointed out in reference 5 that Maxwell’s distribution 
law allows the assumption to be made that rotational 
and translational energies are completely independent, 
so that a molecule with a high translational energy is 
just as likely to have a low as a high rotational energy. 
In the case of vibrational energies, it is shown that those 
molecules that have higher translational velocities also 
carry higher vibrational energies. However, for rea- 
sons that will be presently discussed, the change in 
the vibrational energy is neglected in the energy 
balance equation; hence, we need not be concerned 
with its magnitude in this somewhat simplified analy- 
sis. 


Molecular Energy Transport from the Body 


To determine the energy carried from the body by 
molecules returning to the main mass of gas, it js 
necessary to consider the mechanism of energy transfer 
at the surface of the body. A complete discussion jg 
beyond the scope of this paper; the reader is referred 
to reference 5 for a detailed description of the phe- 
nomena involved. Nevertheless, the more important 
concepts are outlined below because an understanding 
of them is necessary in succeeding computations. 

It appears that, in the practical case, two types of re- 
flection occur. The first and predominating type is 
diffuse reflection. Several subtypes may be differen. 
tiated; however, the most reasonable assumption is 
that the molecules leave the same element of surface 
that they strike and that, in leaving the surface, they 
are emitted in a perfectly random fashion, both as to 
speed and direction and, hence, have a Maxwellian 
velocity distribution corresponding to some as-yet- 
unspecified temperature. The other type of reflection 
is specular in which the molecules leave the body at the 
same angle at which they impinge upon the body. It 
seems logical that diffuse reflection would predominate 
as, in reality, an impinging molecule collides, not with 
a smooth hard surface but with an individual molecule 
or atom of the body material. It would seem probable 
that collisions of this type would obliterate, insofar as 
directions are concerned, any prior history of motion, 
especially if more than one collision occurred at the 
surface. 

As yet, the temperature at which the molecules leave 
the surface has not been specified. In general, the 
molecules that leave the body will not possess an aver- 
age velocity corresponding to the plate temperature. 
In the case with which we are concerned—namely, that 
of a cooled plate—the molecules will leave at some 
temperature that is higher than the plate surface tem- 
perature. This inequality of temperature (or energy 
levels) is defined by introduction of the accommoda- 
tion coefficient, defined as 
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a = — E,)/(E; — E,) (23) 


If all the velocity distributions concerned are Max- 
wellian, the accommodation coefficients may be inter- 
preted in terms of temperature; thus, 


a = T,) 


The accommodation coefficient may be considered as af; 
term that accounts for the inability of the impingingpj 


molecules to adjust themselves to the plate tempera 
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ture during the time they are in contact with the platefsorption 
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nition of the accommodation coefficient given by Eq. 
93), the effects of specular and diffuse reflection on the 
mergy exchange are included in the value of a. 

The value of the accommodation coefficient has been 
measured for several gases in contact with many differ- 
ent materials. Values reported range from less than 
(1 to 1.0, as shown in references 6 to 11. Apparently, 
n0 measurements have been made of the accommoda- 
tion coefficient in a high-speed gas stream. It would 
appear certain, however, that the accommodation co- 
dicient would be directly connected with the average 
time that a molecule spends on the surface—the so- 
ulled time of lingering. At any given instant, there 
must be a certain number of molecules adsorbed on the 
urface, Since it would be impossible to exchange energy 
in zero time. If m molecules per second per unit area 
trike the plate, the number of adsorbed molecules on 
the surface, v is given by v nr where 7 is the time 
interval of adsorption. If the number of molecules on 
the wall is constant or decreases with the mass velocity, 
iswould seem likely, then 7 would decrease as the mass 
velocity increases, and, hence, the accommodation co- 
diicient would be expected to decrease with speed. 

The accommodation coefiicient, if defined rigorously, 
yould have different values for translational, rotational, 
and vibrational energies. However, a series of experi- 
nents performed by Knudsen and described in refer- 
ace 5 showed that the accommodation coefficients 
were numerically identical for rotational and trans- 
ational energies. The vibrational components of the 
ternal energy require much more time to adjust to 


al, oi values than do the rotational and translational 
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wmponents. A discussion of the time of adjustment 
if vibrational energy and numerical values for several 
gases may be found in reference 12. In succeeding cal- 
tulations, it will be assumed that the vibrational com- 
nent of the internal energy is unchanged during the 
hiort time interval that a molecule is in contact with 
he body. By the same line of reasoning, dissociation 
i polyatomic molecules is neglected, since dissocia- 
fon is a function of the intensity of the molecular vi- 
rations. These assumptions are serious in that vi- 
= energy and the energy absorbed during dis- 
sciation would, if present, have a marked effect on the 
te temperature. If the length of time a molecule is 
he on the plate is equal to, or greater than, the 
me required for fully developed vibration to be ex- 
ited, then dissociation would probably occur at the 
figh plate temperatures calculated for high-speed flight. 


hsorption times in a flowing gas. 
To calculate the energy abstracted from the plate due 
)the return of molecules to the stream of gas, use is 


s. Ap-ade of Eq. (23). 


s of the 


>ctional 


E, = E,(1 — a) + aE, (25) 
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iE, is the energy of the molecules if they were at 
bate temperature, then Eq. (25) may be written for a 


h more complete analysis awaits experimental data on 
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monatomic gas that obeys the perfect gas laws as 


E, = E, (1 — a) + onkT, (26) 


and for a diatomic gas that obeys the perfect gas laws 
neglecting vibrational energy as 


= (E, + nkT ;) (1 a) + 5/5 ankT,, (27) 


Radiant Energy Exchange with the Body 


Any body is constantly receiving and emitting radi- 
ant energy. In the steady state, a balance is soon 
reached where the emitted energy is equal to the inci- 
dent energy; the body temperature adjusts itself to a 
value so that this equilibrium is maintained. 

In the present case, the body emits radiation at a rate 
given by 


R, = eoT,* (28) 


The body may absorb radiation by several processes. 
Probably the most important is that caused by direct 
solar radiation. The mean value of the solar radiation 
at the top of the earth’s atmosphere, the solar constant 
Jo, is shown in reference 13 to have a maximum value 
of 93.4 ft.Ibs./sec.-ft.? for a plate normal to the incident 
radiation. The body may also absorb radiation from 
constituents of the atmosphere. Most of the atmos- 
pheric gases do not emit appreciable radiation; how- 
ever, water vapor, ozone, and carbon dioxide emit 
strongly at certain wave lengths. In the interest of 
simplicity, radiation from the atmospheric constituents 
will be neglected, and, where incident radiation is con- 
sidered, it will be assumed to be due solely to solar 
radiation, having its previously mentioned maximum 
value of 93.4 B.t.u./sec.-ft.2.. The amount of this 
energy that is absorbed is, then, 


R, = (29) 


Energy Exchanges Occurring Within the Body 


It is conceivable that, in a missile designed for sus- 
tained flight at very high speeds, a skin cooling system 
could be used to maintain structure temperatures within 
specified limits. Disregarding a discussion of the pos- 
sible cooling systems, we may simply assume that the 
cooling system removes heat at a rate Q ft.lbs./sec.-ft.’. 


Method of Calculation 

Having outlined the processes by which energy is ex- 
changed with the body, it is now possible to calculate 
the temperature of a plate moving at any speed through 
a rarified atmosphere. We may substitute Eqs. (26) or 
(27), (28), and (29) in Eq. (1) and obtain, for a mon- 
atomic gas, where E; = E,, 

aE, — ankT, + — Q=0 (30) 

and for a diatomic gas, where E; = E, + nE,, 


a(E, + nkT) — ankT, + €Jy — eo T,4 =0 
(31) 
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In Eqs. (30) and (31), the value of Z, may be taken from 
Eq. (15) and the value of m from Eq. (6) for the front 
side of the plate. For the rear side of the plate, Eqs. 
(18) and (8) are substituted for values of E,’ and n’, 
respectively. 

For a mixture of monatomic and diatomic gases, the 
incident energies of the components of the gas are 
added together, and the re-emitted energies of the com- 
ponents of the gas are added together, so that the 
energy equation assumes the form: 


LE; — ZE, + eJo — eT, ha Q = 0 


To determine the temperature of a plate with no cool- 
ing except that due to radiation, the factor Q in Eqs. 
(30) and (31) is set equal to zero, and the equations are 
solved for values of 7,. 

To determine the effect of internal cooling on the 
plate temperature, arbitrary values of Q are substi- 
tuted into Eqs. (30) and (31), and the equations are 
again solved for values of 7). 


Scope of Calculations 


Calculations were made to determine the tempera- 
ture of a plate traveling at velocities of from zero to 
20,000 ft. per sec., at altitudes of from 75 to 150 miles. 
The calculations were restricted to plate angles of from 
0° to 20°. The plate, in the first case, was assumed to 
be insulated between the front and rear surfaces so 
that the two surfaces could be treated independently. 
In the second case, the front and rear sides were as- 
sumed to be in perfect thermal contact—that is, the 
temperature of the front side of the plate was equal to 
that of the rear side of the plate. In these two sets 
of calculations, it was assumed that solar radiation was 
absent. 

The above calculations were then repeated, for one 
altitude, flight speed, and plate angle, with arbitrarily 
assumed values of internal cooling and with solar radi- 
ation absent. 

A fourth set of calculations was made for two alti- 
tudes at one value of flight speed and plate angle, con- 
sidering that the plate was uncooled and that solar 
radiation was present. 

The atmospheric properties used in the calculation 
were taken from references 2 and 14. Reference 14 was 
used to obtain the properties of the atmosphere at 75 
miles, while reference 2 was used for the higher alti- 
tudes. There is some difference in the assumed atmos- 
pheric compositions between the tables of references 
2 and 14. Warfield’ assumed that the composition 
by weight was constant, at its sea-level value with in- 
creasing altitude, except for dissociation of oxygen. 
The amount of dissociation was assumed to increase 
with altitude, complete dissociation occurring above 
altitudes of 75 miles. On the other hand, Johnson and 
Slack assumed sea-level composition at all altitudes, 
neglecting dissociation entirely. In order to compare 


JOURNAL OF THE AERONAUTICAL SCIENCES—JULY, 


1948 


values of skin temperature at the several altitudes for 
which computations were made, it was arbitrarily de- 
cided to use the values of temperature and density given 
by Johnson and Slack and to assume complete dis- 
sociation of the assumed constant sea-level percentage 
of oxygen. 


DISCUSSION 


In the computation of the incident energy due to 
molecular translation, it was shown that both the 
number of molecules striking unit urea of the body and 
the resultant kinetic energy incident upon the body were 
functions of the term (U/V,,) sin 6. An examination 
of the graphs of the dimensionless groups x and x’ y, 
and y’, which govern the number of molecules striking 
a section of the body, and the magnitude of the inci- 
dent kinetic energy reveal several interesting features. 
With respect to the number of molecules that strike a 
body on the front side, it may be seen from Fig. 3 that 
this number increases practically linearly with flight 
speed for a given plate or body angle. On the other 
hand, as shown in Fig. 4, the group, x’, that controls 
the number of molecules that strike the rear side of a 
body decreases rapidly with flight speed and is prac- 
tically zero for values of (U/V,,) sin 6 greater than 2. 

From Eqs. (15) and (18) it may be seen that the 
total incident translatory kinetic energy term is com- 
posed of two terms: the first due to mass motion of the 
gas relative to the body and the second due to the 
translatory molecular motions. In the case of the 
front side of the body, the molecular energy term may be 
seen, from Fig: 5, to increase rapidly at first with in- 
creasing values of (U/V,,) sin @ and then to approach a 
constant value of °/, at values of (U/V,,) sin @ greater 
than 2. On the rear side of the same body, the molecu- 
lar energy term continues to decrease rapidly as 
(U/V ) sin increases. 

A simple physical interpretation of the trends noted 
above may be given. First, it is necessary to understand 
that, with a Maxwellian type of velocity distribution, 
most molecules have speeds close to the average molecu- 
lar speed. An extremely small fraction have speeds 
much greater than the average or much less than the 
average. For example, only about 2 per cent of the 
molecules in a given volume of gas have speeds that 
exceed twice the average molecular speed. Thus, it 
should be apparent that at high resultant speeds (the 
resultant speed being the component normal to the 
body, U sin 6) the relatively small spread of molecular 
speeds predicted by the Maxwellian distribution be- 
comes unimportant in comparison with the resultant 
speed. 

This diminishing importance of the spread of molecu- 
lar velocities is the reason for the constancy of the 
second term of Eq. (15) for values of (U/V,,) sin é 
greater than 2. It is interesting to note that the kinetic 
energy term in Eq. (15), which accounts for the trans 
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latory molecular motions, has values ranging from 
2kT at zero velocity to */, kT at values of (U/V,) sin 6 
of 2 and higher. This may be compared with the con- 
stant value of */2 kT given in the simplified analysis of 
Sanger. The larger value at zero velocity is caused by 
the fact that the faster moving molecules are both 
more likely to collide with the body and also carry 
more energy than do the slower speed molecules. The 
actual difference in the total kinetic energy term com- 
puted by means of the simplified theory and the more 
exact theory depends upon the relative magnitudes of 
the terms mU?/2 and WkT. For most of the cases 
computed in this report, the term mU?/2 became large 
in comparison with the term, YkT, at flight speeds 
above 8,000 ft. per sec.; consequently, little difference 
between the simplified and more exact theory would be 
expected at flight speeds greater than 8,000 ft per sec. 

Physically, it is clear why both the number of mol- 
ecules striking the rear side of the body and the transla- 
tory energy incident in the rear side of the body de- 
crease with increasing values of (U/V,,) sin 6. The 
front of the body acts as a shield and prevents collisions 
of the rear side except with those molecules that have 
absolute velocity components in the direction of flight 
and with magnitude equal to the flight speed or higher. 
Consequently, as the flight speed is increased, we should 
expect that the incident energy on the rear of the plate 
would eventually approach zero. 

Fig. 7 shows the magnitude of plate temperatures 
calculated for an altitude of 75 miles in the absence of 
solar radiation. The marked effect of emitted radiation 
from the body in lowering the plate temperature is 
striking. For a large range of flight speeds and small 
plate angles, the plate temperatures are actually lower 
than the ambient air temperature. The curves of tem- 
perature of the rear side of the body show an eventual 
decrease with flight speed in accordance with the dis- 
cussion of the preceding paragraph. It would appear 
that the cooling problem could be alleviated by en- 
suring thermal contact between portions of the body 
inclined at positive and negative angles with respect to 
the flight path. Practically, this could be accom- 
plished by ‘“‘boat-tailing’’ the body and using a heavy 
skin of high thermal conductivity between the nose 
section and the boat-tailed section. It is also apparent 
that, from a skin temperature standpoint, large posi- 
tive body angles should be avoided; in other words, 
the body should have a high fineness ratio. Fortu- 
nately, this dictate coincides with aerodynamic require- 
ments. The effect upon body temperature of perfect 
thermal contact between the front and rear sides of a 
flat plate is shown in Fig. 8. It will be noted that the 
plate temperatures are considerably lower than those 
calculated for the front side of an insulated plate, as 
shown in Fig. 7. 

The effect of altitude on the temperature of the front 
side of an insulated plate is shown in Fig. 9. It is ap- 
parent that the problem of high skin temperatures be- 
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Fic. 8. The effect of flight velocity and angle of attack on 
the temperature of an uncooled plate with front and rear sides in 
perfect thermal contact. 
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comes of negligible importance above altitudes of 100 
miles for flight speeds up to 20,000 ft. per sec. and with 
solar radiation absent. 


|_| 
@=20° 
PLATE | 
Tp, “Fabs AMBIENT AIR TEMPERATURE-z 
| | 
| 
| | | 
200 4 4 
| 
| 
1000 + > + + 4 x 
PLATE | 
TEMP. | | 
| 


390 JOURNAL OF THE AERONAUTICAL SCIENCES—JULY, 1948 


FLIGHT VELOCITY, U=20,000 ft/sec. 
PLATE ANGLE, @=10° 

ACCOMODATION COEFF. = 1.0 
ALTITUDE = 75 MILES 


1200 REAR SIDE INSULATED 


800 - 

PLATE AMBIENT TEMPERATURE 
TEMP., 
Tp, “Fabs. 
| 

400}————_ : 
| 


Q, Btu/sec, sq ft 


Fic. 10. The effect of internal cooling on the temperature of the 
front side of a plate. 
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Fic. 11. The effect of accommodation coefficient on the tem- 
perature of the front side of a plate. 


2000 RADIATION PRESENT 
—SOLAR RADIATION ABSENT 
—- LIGHT VELOCITY, v= 20,000 ft/sec 
| PLATE ANGLE, @ = 10° 
ALTITUDE 
PLATE } 
TEMP, 4 
Tp. *Fobs | 
800 4 
ILES 
| | 
03.04 06 07 08 09 10 
EMMISSIVITY €, DI 


MENSIONLESS 


Fic. 12. The effect of emissivity on the temperature of the 
front side of a plate. 


The effect upon body temperature of internal cooling 
is shown in Fig. 10 for a flight speed of 20,000 ft. per 
sec. and a body angle of 10°. The body temperatures 


are seen to decrease rapidly with the amount of internal 
cooling employed. 

The effect of the magnitude of the accommodation 
coefficient upon the plate temperature at one speed, 
altitude, and plate angle is shown in Fig. 11. It may 
be noted that as the accommodation coefficient is 
lowered, implying increasing difficulty of molecular 
energy exchange, the plate temperature likewise de- 
creases. This phenomena may offer a method of cool- 
ing, since it may be possible to decrease the accommo- 
dation coefficient by altering the surface character- 
istics of the skin. 

The variation of plate temperature with the emis- 
sivity of the plate material in the absence and in the 
presence of solar radiation is shown in Fig. 12. A not- 
able difference between the conditions at 75 and 150 
miles altitude is apparent. At 75 miles altitude, the 
presence of solar radiation has little effect upon plate 
temperatures, the implication being that most of the 
heat transfer is due to molecular energy exchange. The 
plate temperatures rise as the emissivity is lowered be- 
cause of the increasing inability of the body to radiate 
energy arising from molecular collisions. At 150 miles 
altitude, however, the presence of solar radiation 
greatly increases the plate temperature, and, with solar 
radiation present, the plate temperature is practically 
independent of the emissivity. This surprising trend 
may be explained as follows: At 150 miles altitude, the 
magnitude of the molecular energy interchange is small 
in comparison with the radiant energy exchange. Thus, 
neglecting all energy except radiant energy, an energy 
balance on the plate may be written as eJ) = eo7;,’. 
The emissivity € is seen to cancel out, and the plate 
temperature is thus dependent only on the magnitude 
of the solar radiation Jo. Computations of plate tem- 
peratures have not been made with solar radiation 
present for all speeds and altitudes; however, it is 
apparent that at the higher altitudes the magnitude of 
solar radiation becomes increasingly important and the 
magnitude of the emissivity becomes of less importance. 


CONCLUSIONS 


The following conclusions may be drawn from the 
trends shown in the preceding analysis. 

(1) With solar radiation absent, the skin cooling 
problem associated with high-speed missile flight be- 
comes negligible at altitudes of 100 miles and higher up 
to steady flight speeds of 20,000 ft. per sec. and with 
body angles of less than 20°. 

(2) The effect of solar radiation on skin temperatures 
is small at 75 miles altitude at a flight velocity of 20,000 
ft. per sec. but becomes increasingly important as the 
flight altitude is increased. At an altitude of 150 miles, 
solar radiation is the predominating factor that deter- 
mines skin temperature. 

(3) With solar radiation present, the effect of emis- 
sivity on body temperatures becomes of decreasing im- 
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portance as the altitude is increased. At an altitude of 
150 miles, the body temperature, at a given flight speed, 
is independent of the emissivity and depends only on 
the intensity of the solar radiation. 

(4) In order to minimize skin temperatures, positive 
body angles should be kept as small as possible. This 
may be done by designing the body to have high fine- 
ness ratio. 

(5) Skin temperatures may be reduced by ensuring 
thermal contact between portions of the body inclined 
at positive and negative angles with respect to the 
flight path. As much surface as possible should be in- 
clined at negative angles. Practically, this may be ac- 
complished by boat-tailing the body. 
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The Testing of Rotors for Fatigue Life 


JONATHAN WINSON* 
Prewitt Aircraft Company 


SUMMARY 


The object of this paper is to describe an experimental method 
for the fatigue testing of articulated rotor blades. The method 
consists in the application of specified second harmonic control to 
a rotor revolving on a stationary whirl stand. It is shown that 
such a whirl test can induce blade fatigue stresses that well ap- 
proximate those met in flight. Equations are derived to specify 
the values of the test parameters that will allow the duplication of 
a general set of flight fatigue stresses. The test specification for 
a sample rotor is calculated. Some of the theoretical and experi- 
mental considerations involved in the method are discussed. 


INTRODUCTION 


tke, FATIGUE TESTING of articulated rotor blades is 
a subject of much current interest in the rotating 
wing industry. A standard method of testing is re- 
quired which will duplicate flight stresses with a good 
degree of accuracy and yet be reasonable in its practical 
demands. The method must be of sufficient accuracy 
such that it can be relied upon as a measure of safety. 
It must be simple enough and involve a short enough 
period of testing time such that the various companies 
and agencies will feel it desirable to test each new de- 
sign before flight. 

This paper investigates the possibility of satisfying 
these conditions with a whirl stand test of the rotor. 
The proposed test may be described briefly. The rotor 
is run at specified values of rotational speed, collective 
pitch, and second harmonic control. A blade damper 
of specified characteristics is used at the drag hinge. 
The test is run for sufficient time to allow the desired 
number of stress reversals (the reversals occur twice 
per revolution). When the rotor has completed its 
run, it may be considered tested for fatigue. The per- 
fection of this method, then, is the goal that is visualized. 

The present investigation restricts itself essentially 
to a solution of the equations of motion of an articu- 
lated rotor blade under second harmonic control, a 
calculation of the stresses induced in this condition, and 
a comparison to flight stresses as found by methods in 
current use. The conclusion is drawn that there are a 
sufficient number of independent variables to allow the 
possibility of a whirl test duplicating the important 
flight stresses. Equations are written for the various 
test parameters, and a theoretical test specification is 
derived in general terms. The relationships between 
parameters and induced stresses are believed to be a 
good first approximation to the truth. It is indicated 
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that experimental work will be necessary to arrive at 
a final test specification. 

The work is divided into four sections: 

(1) Blade Loading in Flight—The various loads and 
stresses induced in flight are considered. Formulas are 
given for the calculation of the major loads, and the 
specific ‘stresses that the whirl tests must reproduce are 
indicated. 

(II) Blade Loading in Whirl Test. (a) The plane of 
flapping.—-Equations are derived for the flapping mo- 
tion of a blade under second harmonic control and for 
the resulting bending moment distribution in the plane 
of flapping. Independence of those variables respon- 
sible for steady and fatigue stresses is proved. Formu- 
las are derived for the collective and harmonic pitch 
settings necessary to approximate the flight distribu- 
tion of bending moment. 

(b) The plane of rotation.—The motion of the blade 
about the drag hinge, and the in-plane bending moment 
distribution are found. A purely viscous damper is 
assumed to act at the drag hinge. 

(IIT) Sample Calculation—A sample rotor is an- 
alyzed. The complete specification for its fatigue test 
is calculated, 

(IV) Discussion.—Three subjects are considered. 
They are: (1) the question of the induced velocity 
field of the rotor and the method of centrifugal relief; 
(2) the possibility of excessive torsional response to 
second harmonic control; and (3) the path of future 
work in the method. 

All aerodynamic derivations are based on the funda- 
mental work of Glauert! and Goldstein.” In this paper 
only the simplest of articulated rotor blades is con- 
sidered—that is, the blade with constant radial dis- 
tribution of aerodynamic and elastic characteristics. 
The blade is further considered to have the elastic 
center, the center of mass, and the center of pressure 
coincident at each blade section. There is no theo- 
retical barrier, however, to the extension of the analysis 
to more varied designs. 


NOMENCLATURE 

R = radius of the blade 

r = coordinate distance along the blade span 

x = 7/R 

y = coordinate distance perpendicular to blade 
span, defining position of deflected blade with 
respect to its unbent position 

c = chord length 

y = azimuth position of the blade measured in the 


direction of rotor rotation. The blade in its 
downwind position is at zero angle of y 
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Fic. 1. Bending moments and forces at a blade section. 

v = velocity of flow induced through the rotor disc 

V; = component of the resultant velocity of the air 
relative to the blade in a plane perpendicular 
to the blade span 

Vry = component of V; perpendicular to the blade 
span and lying in the plane formed by the 
blade span and the rotor axis of rotation 

Vrz = component of V- which is perpendicular to both 
V,, and the blade span 

dT = thrust component of resultant aerodynamic 
force on segment of blade drdT is parallel to 
Very 

dL = element of lift 

dD = element of profile drag 

a = absolute angle of attack of the blade element 

@ = the induced angle; the angle between V; and Vyz 

w = speed of angular rotation of the rotor 

p = density of the atmosphere 

Qe = slope of the lift curve for infinite aspect ratio 

Cr = coefficient of lift of the blade element 

6 = geometrical blade angle, a function of y 

8 = component of @ independent of y 

6’ = component of @, which is a second harmonic 
function of y 

8 = coning angle 

= lag angle 

= dB/dp,dy/dy 

8,7 = dB/dt, dy/dt 

P = damper constant 

p’ = line density of blade 

E = modulus of elasticity 

I = moment of inertia 

Az, Bz, Cz = spanwise distribution functions 

r = circulation 

B = number of blades 

= tip speed ratio 

M = bending moment 

e = drag hinge offset 

6 = profile drag coefficient 

Subscripts 

T = thrust 

I = inertia 

CF = centrifugal force 

D = drag 

d = damper 

W = weight 

CR = Coriolus 
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(I) LoaDING IN FLIGHT 


Blade stresses m flight are caused, in the last analy- 
sis, by complex distributions of air, inertia, and cen- 
trifugal forces. The resulting load system may be 
schematically portrayed by three components of bend- 
ing moment and three components of force acting at 
each section of the blade. In Fig. 1, a section of blade 
is shown under this general loading. A right-hand 
convesition is used for the moments. 

In common terminology: 


M, = bending moment in the plane of flapping 
M, = torsional moment 

M; = in-plane bending moment 

F, = in-plane shear 

F, = centrifugal tension 

F,; = shear in the plane of flapping 


The condition of interest in the fatigue analysis of a 
rotor blade is that of forward flight at a velocity equal 
to or above cruising speed. In such a flight condition, 
the moments and forces assume characteristic span- 
wise distributions and combine to cause high values 
of fatigue stress at two blade sections. At about the 
mid-span of the blade, the steady centrifugal] stress com- 
bines with the high fluctuating stress due to moment 
M, to give the design fatigue condition in the plane of 
flapping. The fluctuating moment M; contributes, 
secondarily, to the fatigue system at this point. At 
the inboard end of the blade, the fluctuating stress due 
to M; reaches a maximum and combines with the cen- 
trifugal stress to give the design condition in the plane 
of rotation. In the balanced type of blade considered, 
the stresses due to We, Fi, and F2 are small. This 
is the picture given by current methods of rotor 
analysis. 4, 43, and F; will now be treated more 
completely. 


M,, Bending Moment in the Plane of Flapping 


The blade is hinged at the inboard end and may move 
without restraint about the hinge. The characteristics 
of the spanwise bending moment distribution are largely 
determined by the four boundary conditions of the sys- 
tem: 


(1) Deflection of the blade at the hinge = 0. 
(2) Bending moment at the hinge = 0. 

(3) Bending moment at the tip = 0. 

(4) Shear at the tip = 0. 


In the Appendix, equations are given for the quanti- 
tative determination of .4 and the associated shear 
F;.3:4 Specifically, the equations allow the calculation 
of flapping motion in flight and the bending moment 
distribution on an inflexible (but coned) blade. They 
are carried through second harmonic terms. Provision 
is made for the inclusion of first harmonic control and 
an unsymmetrical fore-and-aft induced velocity dis- 
tribution. The question of blade flexibility is treated 
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in section Ila, and formulas are given there for the cal- 
culation of centrifugal relief. 

Considering the inflexible blade bending moment 
distribution at any one azimuth position, the following 
may be seen: All moments are made up of various 
aerodynamic and inertia constants multiplied by the 
three basic spanwise distribution factors 

A, = (1/2) — x + (x?/2) 

B, = (1/3) — (x/2) + (x*/6) 

C, = (1/4) — (x%/3) + (x*/12) 
The aerodynamic moment consists of all three distribu- 
tions; the inertia moments are of the B, distribution 
exclusively. Considering that the resultant bending 
moment is essentially a balance between aerodynamic 
and inertia moments (MM; — M,; — Mer) and, further, 
that the moment at the inboard hinge is zero, it may be 
deduced that the moment is made up of (C, —,B,) and 
(A, — B,) distributions, each with zero value at the 
hinge. - The relative importance of the two distributions 
is proportional to the relative magnitudes of the con- 
stants multiplying C, and A, in the expression for aero- 
dynamic moment. Inspection of the equation for this 
quantity shows that the A, term is of second order by 
the ratio of approximately yp’ to 1. 

The resulting bending moment on the inflexible blade 
is then a (C, — B,) distribution with a second order 
(A, — B,) distribution, both adjusted to zero at the 
blade.root. The (C, — B,) distribution has a peak at 
x = 0.422.’ Fig. 2 and Table 1 define (C, — B,), 
normalized at its maximum value. The (A, — B,) dis- 
tribution reaches its peak at x = '/3. Many calcula- 
tions have established the approximate limiting rela- 
tionship 


(A, B,)/(C, B,) = 0.15 


for normal flight. (4, — B,) may add to or subtract 
from (C, — B,). Using this relationship to gain some 
idea of the variation of spanwise position of the maxi- 
mum value of moment distribution as the blade pro- 
ceeds around the azimuth, it is found that the peak 
moves between x = 0.405 and x = 0.448. The total 
distributions in the limiting cases vary only a small 
amount from (C, — B,). The application of centrif- 
ugal relief causes all peaks to move outboard, the 
amount of movement and final distribution depending 
on centrifugal and elastic constants. The spanwise 
limits of maximum value remain in a range of approxi- 
mately 0.14 blade radius. 

As the blade sweeps around the azimuth, the peak 
value of bending moment varies from minimum to 
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° Fic. 2. Norimalized (C, — B,) distribution. 
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Fic. 3. Variation of maximum moment with azimuth. 


maximum. The variation is first harmonic with an 
appreciable second harmonic contribution. Fig. 3 
gives a typical curve of variation of peak moment with 
azimuth.. 


M:;, In-Plane Bending Moment 


There is a vertical pin at the inboard end of the blade 
about which the blade may oscillate in its plane of rota- 
tion. Restraint is provided by a damper at the pin. 
The boundary conditions of the system are: 


(1) Bending moment at root = moment exerted by | 


damper. 

(2) Bending moment at tip = 0. 

(3) Shear at tip = 0. 

With dampers currently used, the bending moment 
distribution varies rather uniformly from the damper 
moment at the blade root to zero at the blade tip (ac- 
cording to a combination of A,, B,, and C, functions). 
When the variation of moment with azimuth is con- 
sidered, it must be stated that analytical work concern- 
ing the moment is by no means as complete as that 
dealing with the moment 1/,. The damper characteris- 


TABLE 1 
(C, — B,) Distribution 


0.2 0.3 0.4 0.422 
0.690 0.905 0.997 1.00 


M, 0 


0.5 0.6 0.7 0.8 0.9 1.0 
0.962 0.813 0.583 0.318 0.097 0 
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tic usually makes strict analysis difficult. In addition, 
the inboard end of the blade is frequently designed by 
other conditions, and .a careful fatigue analysis is 
avoided by the use of broad and conservative assump- 
tions as to the situation in flight. When an analysis is 
carried out, it is usually assumed that the damper is 
purely viscous in nature, and the lag motion is approxi- 
mated by a first harmonic function. This leads to the 
picture, at the root, of a steady centrifugal stress and a 
harmonic fluctuating stress varying from positive to 
negative maximum once per revolution. The true 
motion undoubtedly has appreciable higher harmonic 
components with the consequent asymmetry of the 
fluctuating moments about zero, 


F,, Centrifugal Tension 


The centrifugal tension gives rise to the major por- 
tion of the steady stress in the fatigue system. It is 
given by: 


Fy, = w*R*p’ x dx (1) 


From this descripiion it may be concluded that a 
successful whirl test must produce, as a minimum, the 
following system of stresses: 

(1) At a point about midspan, a given tension stress, 
a given steady bending stress, and a given fluctuating 
bending stress (both in the plane of flapping). 

(2) At the inboard end, a given steady tensiom stress 
and a given fluctuating bending stress in the plane of 
rotation (the fluctuating stress varying from a positive 
maximum to an equal negative maximum). 

In the next section, it will be shown that these re- 
quirements can be fulfilled, to a good degree of approxi- 
mation, by one whirl test. 


(IIA) Wuirt Test, PLANE OF FLAPPING 


Inflexible Blade 


Fig. 4 is a diagram of the velocities and forces at a 
blade element. With the common blade element as- 
sumptions: 


dL 


Velocities and forces at a blade element. 


Fic. 4. 


dT = (p/2)C,V,.°c dr (2) 
On the test stand: 
(3) 
Also: 
= aoa (4) 
a = 0 — (5) 
Vy =v+rB (6) 


Assuming partially constant and partially triangular in- 
flow distribution : 


v=u+u'(r/R) (7) 


The geometrical blade angle under second harmonic 
control may be expressed as? 


6 = + sin 2y (8) 

Let there be defined : 
v = u/wR (9) 
v’ = u’/wR (10) 


Making the sequence of substitutions of Eqs. (35) 
through (10) in Eq. (2) and recognizing the rela- 
tion B/w = 8’, there is obtained for the thrust distribu- 
tion: 
dT 
wre + 6’ sin — v’ — B’)—r(vR)] 
r 

(11) 


A double integration of the thrust distribution car- 
tied out with respect to span yields the expression for 
the bending moment .distribution due to aerodynamic 
thrust. Thus: 


Mr = (dT /dr) dr dr (12) 
Mr = K [C,(0 + 6’ sin 2y — v’ — 8’) — By] (13) 
where 
K = (p, 2)a .w*cR* 


There are three other sources of transverse bending 
moment—namely, flapping inertia, the transverse 
component of centrifugal force, and blade weight. The 
last is small and is not considered. The bending mo- 
ment distributions due to inertia and centrifugal force 
are now found. 


M, = dr dr = —w*R*p'B"B, (14) 
Mer = dr dr = —w*R*p'BB, (15) 


The resultant bending moment in the plane of flap- 
ping on the inflexible blade is the sum of the three 
component distributions. Thus:. 


M, = C,K(@ + 6’ sin 2y — v’ — 
B,[Kv + w?R*p'(8 + 8”)] (16) 
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The condition of zero bending moment at the root 
applied to Eq. (16) yields the differential equation for 
blade flapping motion. The equation is: 


| 
O = (+ sin — —p’) 


| 
3 [Kv +w?R%p’(8 + B”)] (17) 


The exact steady state solution to Eq. (17) is: 


B = Bo + Bi cos 2y + Be sin 2y (18) 
where 
Bo = (3/4)1[0 — v’ — (4/3)r] (19) 
—6'[P/2(4 + P)] (20) 
= —6'[I/(4 + P)] (21) 
and 
l= K/w*R'p’ (22) 


Eqs. (18) through (21) are substituted in Eq. (16). 
There result two parts to the moment distribution 
M,—namely, Mic, independent of azimuth angle, and 
My, a second harmonic function of azimuth angle. 
The expressions for the two quantities are: 


Mic = K(@% = v’) (3 4)B,]) (23) 


Miy = sin — 3B,) + 


Kl 
cos 2y (6 +P 
Both constant and variable distributions are of the 
form (C, — B,) with zero value at the root. (C, — B,) 
: is the distribution given in Fig. 2 and Table 1. ,It 
“i has its maximum at x = 0.422. Mic and My are 
rewritten in terms of their maximum values. Thus, 
Mic is constant with.azimuth and has the maximum 
value 


) 2C, — (3/2) Be) (24) 


Mic (maximum) = 0.01084 K(6> — v’) (25) 


M,y is a second harmonic function of azimuth and has 
the maximum value 


M,,(maximum) = 0.02168 K#’/V4+F (26) 
which occurs at the position of azimuth defined by: 
tan 2y = 2/1 (27) 


It is seen, therefore, that, as far as the inflexible blade 
bending moment in the plane of flapping is concerned, 
the criterion for duplication of flight stresses is fulfilled. 
There is a constant bending moment distribution de- 
pendent only on  (v’ is shown below to be a function 
of %). Both distributions are of the (C, — B,) type, 
which is very nearly that met in flight. It will be found 
later that the consideration of blade flexibility does not 
change these characteristics. 
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Fic. 5. Approximation to inflow distribution. 


Determination of the Induced Velocity 


Goldstein's Vortex Theory’ gives, for the induced 
velocity field, 


v = BI sin dk (28) 
where k is an inflow factor. Since 
pV PCr (¢/2) = oV,T (29) 
then 
v = BcV,C,/8zxr sin 6k (30) 


Making the common blade element substitutions 


= a, [4 (v/wr)] 
sin = v/wr 


there results the quadratic equation for induced ve- 
locity: 


(31) 


where 
S = 


This theory does not take into account the effect of 
the azimuthal variations of blade angle and coning 
angle. The factor k varies appreciably from 1 only at 
the blade tip. It allows the accurate evaluation of the 
tip loss. The tip loss is neglected in this work, con- 
sistent with the fact that the equations used for mo- 
ments in flight do not include it. These points are 
treated more fully in the discussion. Taking k = 1 then 
and substituting Eqs. (7), (9) and (10) in Eq. (31), there 
is obtained: 


(2/S)(v + v’x) = -1 + V1 + 4(6/S)x (32) 


The accuracy of the assumption of partially con- 
stant and partially triangular in-flow distribution may 
be judged from Fig. 5, where the actual flow and the 
approximation are shown at a typical value of 6/5. 
A chart of v and v’ versus 6/S is given in Fig. 6. The 
use of the chart allows Eq. (25) to be evaluated in a 
straightforward manner. 
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Fic. 6. Variation of inflow ratios with collective pitch. 


The Effect of Blade Flexibility 


There remains the problem of correcting Eqs. (25), 


(26), and (27) for the effect of blade flexibility. The 
differential equation describing the system of a rotating 
elastic blade under transverse loading is: 


* dy 
2Ko d*y dS; dS; (33) 
ae 
where 


Ky = p’w*R*/2EI 


and d.S,/dr is the transverse load on the inflexible blade 
due to aerodynamic, inertia, and transverse centrifugal 
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forces. d5S,/dr is made up of constant and second har- 
monic parts. Thus: 


dS;/dr = fi(x) + fe(x) sin + fe(x) cos (34) 
The exact solution of Eq. (33) is 
= yi(x) + yo(x) sin + ys(x) cos 2p (35) 


Substitution of Eq. (35) in Eq. (33) and separation of 
harmonics lead to the following two characteristic 
equations: 

dy, dy, 


Kell — 29) + fle) (36) 


and, as typical for the harmonic functions: 


d*yo dy. _ 
9 


R* 
SKoy2. = (x) (37) 


Owen® has presented an excellent method for the 
solution of this type of differential equation. The 


method has been modified somewhat by the author‘ and * 


may be summarized as follows: 
A bending moment distribution on the inflexible 


blade given by 


(inflexible) = R A, E + + | 
(38) 


vields a resultant bending moment on the actual blade of 


TABLE 2 
Centrifugal Relief, Steady Part 


Ag function, a7 = —0.038094 


x 0 ae | 0.2 0.4 ° 0.5 0.6 0.8 0.9 1.0 
Mi 0.000000 0.081000 0. 128000 0.144000 0.125000 0.096000 0.032000 0.009000 0 
Me 0.012698 0.012623 0.012168 0.009555 0.007542 0.005350 0.001537 0.000401 0 
Mai 0.333333 0.283500 0.234666 0.144000 0.104166 0.069333 0.018666 0.004833 0 
Me 0.000000 0.000150 0.001082 0.006912 0.011979 0.018288 0.033450 0.041675 0.050000 
Az function, a = —0.004365 
x 0 0.1 0.2 0.4 0.5 0.6 0.8 0.9 1.0 
10° X pe 1.4550 1.4548 1.4492 1.3314 1.1615 0.9063 0.2979 0.0800 , 0 
My 0.000000 0.000810 0.005120 0.023040 0.031250 0.034560 0.020480 0.007290 0 
Mai 0.333333 0.283500 0.234666 0.144000 0.104166 0.069333 0.018666 0.004833 0 
10° X us 0 0.000456 - 0.0121 0.2780 0.7069 1.4442 3.9010 5.4832 7.1429 
TABLE 3 
Centrifugal Relief, Second Harmonic Part 
Ag function 
x 0 0.1 0.2 0.4 0.5 0.6 0.8 0.9 1.0 
Wa —1.0 —0.8505 —0.7041 —0.4320 —0.3123 —0.2079 —0.0561 —0.0145 0 
Me —0.03809 —0.03151 —0.02530 —0.01485 —0.01065 —0.00712 —0.00192 —0.00156 0 
Az function 
x 0 0.1 0.2 0.4 0.5 0.6 0.8 0.9 1.0 
Hs —1.0 —0.8505 —0.7041 —0.4320 —0.3123 —0.2079 —0.0561 —0.0145 0 
—0.004315 —0.003657 —0.002943 —0.001632 —0.001100 —0.00067 —0.00017 —0.00004 0 
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My(actual) = (EI/R*)A,m (39) 
and a deflection curve given by 
y = A,(us + dor) (40) 


The tabular values of the yu coefficients for two of the 
A, functions are given in Tables 2 and 3. Table 2 is for 
the solution of Eq. (36) and is used to modify Eq. (25). 
Table 3 relists those coefficients that are changed when 
the solution is sought for Eq. (37). mw: and ps are the 
same for either solution. Table 3 is used to find the 
effect of flexibility on Eq. (26). 

The application of the method consists in approxi- 
mating the value of inflexible blade bending moment 
given by Eq. (25) or Eq. (26) (taken in conjunction with 
the known (C, — B,) type distribution) by the right- 
hand side of Eq. (38). The A» and Az, functions are 
used separately, or the two are used together in any 
linear combination. The specific values of Ay and A: are 
determined in the process of approximation. These 
values are then substituted in Eq. (39), and the actual 
moment distribution is found. The blade deflection 
curve can be calculated from Eq. (40). 

The method may be used in conjunction with the 
equations in the Appendix to find the true bending 
moments in forward flight. Table 2 is then used and a 
good approximation to the exact centrifugal relief effect 
is obtained.° 

In the present analysis, the (C, — B,) distribution 
has been considered for both constant and second har- 
monic moments. Fig. 7 shows the moment distribu- 
tions that result when the inflexible moment distribu- 
tion of Fig. 2 is analyzed at various values of the quan- 
tity w*R'p’/EI. Two centrifugal relief factors are de- 
fined: 

_ Mic (maximum, inflexible) 


F= 


41) 
Mic (maximum, actual) (41) 


M,, (maximum, inflexible) 


F’ = (42) 


My (maximum, actual) 
Fig. 8 is a plot of the two factors versus the quantity 
w*R‘p’/EI. The spanwise position of the peaks of the 
final distributions are also indicated. 
Combining Eq. (41) with Eq. (25) and Eq. (42) with 
Eq. (26), there are obtained the basic equations for the 


determination of two of the test parameters. Thus: 
(% — v’) = FMi-/0.01084K (43) 
6’ = F’My V4 + 2/0.02168K (44) 


where Mig and My represent the maximum values of 
the constant and variable flight bending moments which 
are to be reproduced. 


(IIs) Test, PLANE OF ROTATION 


Fig. 9 shows the important geometrical relationships 
in the plane of rotation. A purely viscous damper acts 
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Fic. 8. Centrifugal relief chart. 


| 


Fic. 9. Geometrical relationships in the plane of rotation. 


at the drag hinge. There are bending moment dis- 
tributions in the plane of rotation due to the following 
forces: (1) aerodynamic drag, (2) lag motion inertia, 
(3) component of centrifugal force, and (4) reversed 
effective force of the Coriolus acceleration. In addition 
there is a moment due to the damper at the blade root. 
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Aerodynamic Drag 


=) 


A standard blade element analysis! leads to the final 
expression : 


Mp 


6 


y?C, — — K sin 2¥(vB, + v’C,)(0’ + 
261) + K cos 2¥(vB, + v'C,)(2B2) (46) 


Inertia 


y is the lag angle measured positively in the direc- 


tion opposite to that of rotor rotation. Then: 
M, = rp'y dr dr = w?R*p'y"B, (47) 
Centrifugal Force 
For small lag angles: 

Mer = dr dr = d w?R*p'yA, (48) 
where d = e/R. 
Coriolus Force 
Mer = dr = 2p'w?R*B'BB, (49) 


Expanding the product 8’8 there is obtained 
= 2p’w?R®B, [sin 2601) + cos + 


sin + Bs) + cos 4¥ (26:82)] (50) 
Damper 
The damper moment is given by the expression: 
Ma = Py’ = Pwy (51) 


where P is the damper constant. 


The bending moment distribution in the plane of . 


rotation is: 
Ms = Mp + Mi + Mer + Mer (52) 
At the blade root 
= —Pwy’ (53) 


Eq. (53) is the differential equation for y. Its exact 


solution is the function: 


Y= yo + yi cos 2y + yo sin + y3 cos 4y + 
yasin 4y (54) 


in which 
_2l 
(55) 
4 d 4 d\? 
(56) 


(57) 

od 16 d\? 

3 = + (48 - 4) + 
(58) 

(59) 


The subsidiary definitions are: 


P’ = P/wR*p’ 
+ 28) 
= 3 0 By 3 4 ( =P1 


4 , 
= Bus + i(2 + ) 


az = — Bi’) 
ay = (4/3) 


The bending moment of interest is that at the root 
given by Eq. (53). Substituting Eqs. (56) through 
(59) in Eq. (53), it is seen that the bending moment at 
the root consists of a second harmonic part and a fourth 
harmonic part. The value of the maximum amplitude 
of the second harmonic part is 


2 
ep 


The maximum occurs at the azimuth angle defined by 


tan = 61 
an (61) 
a2 + 2 ay 
For the fourth harmonic the maximum is 
| 2 2 
a3? + 
My = 4u*R'p'P’ 
At an azimuth position 
6 d 
+ 4P’as 
3 2 3) 
tan 2y = 6: 
as 3 9 a3 


It will be found in most practical cases that the numeri- 
cal value of Eq. (62) is negligible as compared to Eq. 
(60). It will further be found that Eq. (60) is closely 
approximated by 
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(a? + ay’) 
M; = 2wP | (64) 
2 


3 


The criterion for the duplication of flight stresses 
in the plane of rotation is therefore satisfied. The 
choice of the proper damper constant P will produce a 
symmetrical fluctuating bending moment. of any de- 
sired magnitude. 


(III) SAMPLE CALCULATION 


A rotor of the following characteristics is presented 
for test: 


B =3 

R = 17ft. 

E = 30,000,000 Ibs. per sq.in. 

I = 0.160 in.‘ (plane of flapping) 

c = lft. 

w = 82 rad. per sec. in normal flight 


p’ = 0.0901 slugs per ft. 
a.= 59 

6 = 0.010 

e = 1.5ft. 


It is desired to reproduce, at the point x = 0.5, the 
centrifugal tension due to rotation at flight r.p.m. and a 
bending moment in the plane of flapping varying from 
2,900 to 1,100 in.lbs. At the inboard end of the blade 
it is desired to produce the normal centrifugal tension 
and a bending moment in the plane of rotation vary- 
ing from plus to minus 2,000 in.lbs. (the lag motion 
in the flight analysis having been assumed first har- 
monic). The test is run at flight r.p.m. 


The elastic constant is first calculated. 
w*R'p’ 32? X 174 X 0.0901 X 144 
EI 30,000,000 X 0.160 


= 230.5 


From Fig. 8, F = 6.95, F’ = 3.26, and the peak of the 
bending moment curves occur at x = 0.520. 


Various other constants are found: 


K = 5 = 


0.002378 X 5.9 X 322 X 1.0 X 17% 


‘ = 598,000 

598,000 

w*R3p’ 3392? X 177 X 0.0901 

Beaw 3X10X59 
S= = = 0.04 

SrR Sr X17 

d = — = — = 0.0882 


4+2=4+4 1.74 = 5.74 
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The mid-span fatigue system consists of normal cen- 
trifugal force, steady bending moment 2,000 in.lbs., 
and fluctuating moment plus and minus 900 in.lbs, 
The operation of the test at flight r.p.m. repro- 
duces the proper centrifugal tension. It is noted that 
the peak of the test bending moment curve occurs at 
0.520 rather than 0.500. It is decided to reproduce 
accurately the condition at x = 0.500. This means that 
a portion of the blade otitboard of 0.500 will be stressed 
to a slightly higher extent than in flight, while a portion 
inboard will be slightly understressed. In a more care- 
ful application of the method, the designer adjusts this 
small effect to best suit what he desires of the test. 


Then, for the steady bending moment, Eq. (43) gives 
(00 v’) = FM,,-/0.01084K 


for a peak at 0.520. Estimating an adjustment from 
the blade moment values at x = 0.515and 0.495, w?R‘p’ + 
EI = 200, Fig. 7, a factor is arrived at which is so close 
tolthatlisused. Therefore: 


, 6.95 X 2000 
— 


= : = 0.1788 
0.01084 & 598,000 x 12 


Using Fig. 6, it is seen that % = 0.239 (13.7°) satisfies 


this relationship. Subsidiary values are: 
vy’ /S = 1.47, v’ = 0.0608 
v = 0.0207 


'S == 5.78, 
= 0.5, 


Considering the fluctuating moment, Eq. (44) gives: 


x = 0.0453 or 2.6° 
0.02168 X 598,000 X 12 


The flapping coefficients are : 


3 4 
Bo = = 3/4 X 1.32 (0.1788 — 


0.0276) 
Bo = 0.1493 or 8.55? 
1.74 
= —0.0453 (524) = —0.00687 
Bi X 574 
1.3 


For the inboard end of the blade: 


a = (4/3) X 0.1493 & —0.00687 + 
1.32(0.0221) (0.0315) 
a, = —0.00045 
a, = (4/3) X 0.1493 K —0.0104 + 
1.32(0.0221)(— 2) (0.0104) 
a2, = —0.00268 


[(4/3) — (d/2)]? = 1.289? = 1.66 
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TESTING ROTORS FOR FATIGUE LIFE 


From Eq. (64): 
12 X 2 X 32.0 ¥ 0.00000738 
P = 1,235 ft.lbs. sec. 


To summarize, the test is specified by the four quanti- 
ties: 


w = 32 rad. per sec. 
= 
6’ = 2.6° 
P = 1,235 ft.lbs. sec. 


(1V) Discussion 


There are several points that require further analysis. 
The hydrodynamic problem of the velocity field induced 
by a C, distribution variable with azimuth can be at- 
tacked theoretically but is beyond the scope of this 
paper. Concerning the error involved in neglecting 
the cyclic effect, it may be noted that the variation of 
C, with azimuth is not severe. There is usually an am- 
plitude of variation of some 20 per cent of the average 
value. 

The introduction of tip loss effects has been considered 
an excessive refinement for this first analysis. The 
effect of tip loss is appreciable, however, and its 
introduction will be required when the testing method 
is considered more deeply. Whatever theory of tip 
loss is assumed must be consistently applied to both 
flight and test equations. 

Concerning the method of centrifugal relief, although 
the formulas used are of a high degree of refinement, 
aerodynamic damping due to blade flexing is not con- 
sidered. It is expected that the inclusion of this factor 
will have a small effect on the value of F’ and a some- 
what larger effect on the azimuth angle at which the 
maximum moment occurs. 

The question of torsional response to harmonic 
control is of importance. The test, with control at twice 
rotor speed, undoubtedly brings the vibratory system 
closer to resonance than does single harmonic control. 
The lowest natural frequency of the system is usually 
sufficiently high, however, so that this effect remains of 
reasonable proportions. In addition, a tool for the 
elimination of the difficulty exists in the fact that a 
stiff and damped test stand control system can replace 
the flight control installation. 

It should be indicated that there is no unique quality 
associated with second harmonic control as opposed to 
control at other multiples of rotor speed. The use of 
second harmonic control is a compromise. It is, as far 


as can be seen theoretically, the most satisfactory mean 
between control at lower multiples of rotor speed, 
which leads to excessive angles of attack in the at- 
tempt to reproduce the required stress pattern, and 
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control at higher multiples of rotor speed, which ac- 
centuates the danger of torsional vibration. 

To summarize, this has been a first consideration of 
the concept of a whirl stand fatigue test for rotors. 
Much work remains to be done, both theoretical and 
empirical. A precise sequence of experiments is re- 
quired to check the method, modify it, and bring it 
into the practical realm. There is, in addition, the ex- 
tension of the work to teetering and unarticulated 
rotors. It is hoped that the necessary efforts will be 
applied in these directions, for, should the method prove 
successful, an important step will have been taken 
toward the goal of reliability in the rotating wing 
machine. 


Appendix* 


EQUATIONS FOR FLAPPING MOTION 
Let 
B = da — cos — a sin — a3 cos 2W — a, sin 2y 
Define cyclic control as follows: 
Aé = 6 Sin ¥ + 6 cosy 
Further define the following constants: 


Ar = (60/4) (1 + + — (v/3) 


As = p/3 

= —p?/8 

As = (0:/4) (1 + */om?) + — 
As = — — (u?/2)] 

Au = — u/6 


Ais = (62/4) [1 + (u*/2)] — (1/3) 


Au = 1/4[1 + (u*/2)] 
Ais = (6024/3) — (uri /4) 
Au = — (Bou?/ 4) — (0qu/3) 


K, = pa,R'c/2Iy 


Then the flapping coefficients are given by the solu- 
tion of the following five algebraic equations: 


— (ao/Ki) — ayy = Wor / Kyl yw? (la) 
As + aids + adn = 0 (2a) 

— GoAz + + = O (3a) 

Ais — — (43/2) — (3a4/Ki) = (4a) 
— — (3a3/Ki) + (a4/2) = (5a) 


* All equations in this section are reduced from a set dealing 
with tapered and twisted rotor blades with an arbitrary degree 
of flapping-feathering coupling. The original derivations are to 
be found in references 3 and 4. 
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M,, BENDING MOMENT ON AN INFLEXIBLE BLADE 


M, = Mr M, Mer My (6a) 


The equations for these quantities are as follows: 
Mr = KC,[9% B’ + Ay cos y + sin + KuB, = 3) cos y + (20 = B’) sin y 


1 : 1 0. . 
6, cos + sin 2Y + (—3a3 + a4) cos + (—3a, — az) sin 3y | + Ky?A, 2 + cos 


3 8 
(2, sin y — 608 2y sin + cos 3y sin + cos 4y — (7a) 


M, = p’w*R®p"B, (8a) REFERENCES 
‘ Mer = p’wR®BB, (9a) 1 Durand, W. F., Aerodynamic Theory, \st Ed., Vol. 4, Chapt. 
{ Mw = W,RA, (10a) 10; Julius Springer, Berlin, 1934. 


2 Kaman, C. H., Aerodynamic Considerations of Rotors in Hov- 


The additional quantities involved in these equations 
ering and Vertical Climb Conditions, Journal of the Aeronautical 


a icht of the blad Sciences, Vol. 10, No. 2, p. 201, 19438. 
. I oe fj " ; f the blad iis h 3 Winson, J., Equilibrium of a Helicopter in Forward Flight, 
: | he moment 0 ‘merGa 0 e ade about the Kellett Aircraft Company Report No. 230.32, December, 1944. 
a flapping y= F 4 Winson, J., Bending Moments of a Helicopter Blade in Forward 
? = coordinate of the center of gravity of the, Flight, Kellett Aircraft Company Report No. 250.28, July, 1945, 
blade , : 5 Duberg, J. E., and Luecker, A. R., N.A.C.A. Wartime Report 
1 = m/wR where ™ is the maximum value of the No. L5E23, Washington, D.C., August, 1945. 
fore and aft triangular inflow distribution in 6 Owen, J. B. B., The Stressing of Gyroplane Blades in Steady 
forward flight Flight, R. & M. No. 1875, British A.R.C., 1939. 
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SUMMARY 


Test results indicate that the angle of attack becomes critical 
for relative inlet velocities in excess of 1,400 ft. per sec. for air 
and at comparable Mach Numbers for the various Freons and 
that the optimum positive angle of attack is approximately 4°. 

The curves shown in Fig. 9 indicate that, particularly for low 
inlet temperatures, the increase in inlet velocity due to chartge 
of density becomes marked for approach velocities in excess of 
400 ft. per sec. 


INTRODUCTION 


_" PRESENT ARTICLE is intended to facilitate the 
work of the practical designer who is faced with 
the problem of designing a highly efficient centrifugal 
impeller operating with high inlet flow velocity. As 
this in effect is the equivalent of high inlet Mach Num- 
ber, most of the comments are equally applicable to the 
newer designs of industrial centrifugal compressors 
with electric drive and step-up gear that handle Freon 
vapor in a refrigeration cycle, since the velocity of sound 
in F-113 is less than 400 ft. per sec. at some operating 
temperatures. 

While the effect of change of density as a result of 
acceleration may be neglected in preliminary calcula- 
tions for all but extreme cases of high inlet velocity, 
the final design must take note of this effect, and for 
ready reference a step-by-step analysis of the method of 
constructing curves to give this increase in velocity is 
presented. 


INLET DESIGN CONSIDERATIONS 


In designing high-speed centrifugal impellers such as 
are used in turbosupercharger compressors, geared air- 
craft engine superchargers, and gas turbines, the theo- 
retically correct impeller inlet vane angle 8, shown in 
Fig. 1 is easily obtained. 

A typical inlet velocity diagrams is shown in Fig. 1, 
where C, is the true axial approach velocity of the flow, 
corrected for the acceleration due to change in density 
of the air at the impeller inlet, as described later; Uj is 
the tip, or eye, velocity of the maximum inlet diameter 
of the impeller; the closing side of the triangle gives 
the magnitude and direction of the relative inlet ve- 
locity W, of the flow, and theoretically the vane inlet 
angle should coincide with the direction of W,. 
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Inlet Conditions of Centrifugal Compressors 
for Aircraft Engine Superchargers 
and Gas Turbines 


H. WOODHOUSE* 


General Engineering § Research Corporation 


Furthermore, it may be seen from the composite 
velocity diagram of Fig. 2 that, theoretically, the 
tangent of the vane inlet angle should vary inversely 
with the peripheral velocity of the vane inlet edge. 

If U, and U, are the peripheral velocities at two arbi- 
trary inlet diameters (see Fig. 3), then, since C; is a con- 
stant value over the.entire impeller inlet area, we have: 


= U, tan = U, tan 8, 


In practice it is necessary to increase the vane inlet 
angle several degrees over its theoretical value for the 
impeller to develop its maximum efficiency at the de- 
signed flow. In fact, if the inlet angles of some high- 
speed, high-performance impellers are not increased in 
this manner, they will completely fail to pass the re- 
quired flow. 

With all existing high-speed ceutrifugal impeller 
types, the air is made to undergo a sudden change of 
direction immediately after entering the impeller inlet 
(see Fig. 4), which sets up considerable turbulence in the 
flow, causing a “choking” effect and preventing the 


U; tan By => 


- flow from utilizing fully the available inlet area, thereby 


making an increase in angle necessary in order to pro- 
vide an effective area equal to the designed value. 

This increase in inlet angle (see Fig. 5) over its theo- 
retical value is known as the “‘angle of attack.” 

From any series of test runs at different inlet tip 
speeds Uj, the inlet flow velocity C; may be calculated 
for a number of different flows at each speed and the 
corresponding theoretical inlet angle 8; may be ob- 
tained, as in the velocity diagram of Fig. 1. By de- 
ducting this value of 6, from the actual vane angle in 
each case, the angle of attack for any particular speed 
and flow may be obtained. 

If, for any one speed, the values obtained for the angle 
of attack for a wide range of flows are plotted against 
efficiency, as in Fig. 6, it will be seen that maximum 
efficiency occurs.at one particular value of the angle of 
attack. It will also be observed, if similar curves are 
plotted for the complete range of speeds, that at low 
speeds the curves are flat, indicating no particular 
critical value of the angle of attack. However, as the 
speed is increased, the curves become steeper until at 
the highest speeds they are almost vertical, giving only 
one value for the angle of attack. 

This suggests that, as the relative inlet flow velocity 
approaches the supersonic region (or, in other words, 
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INLET CONDITIONS OF CENTRIFUGAL COMPRESSORS 405 


' the higher its Mach Number), the closer must the 


actual inlet vane angle conform to the correct angle of 
attack for efficient operation. Recent tests by the 
author on centrifugal Freon refrigerating compressors 
running at lower tip speeds but comparable Mach 
Numbers indicate that the pattern of the curves is sus- 
tained for centrifugal compressors of conventional in- 
dustrial design. 

Only tests with vaneless diffusers should be used for 
the determination of the angle of attack, since the 
flow might be affected by the characteristics of a vaned 
diffuser. 

The influence of the angle of attack will be much more 
noticeable with impellers of high, than with those of 
low, specific speed (see Fig. 7), since, if both have the 
saine outside diameter tip speed, the high specific speed 
impeller will have a relative approach velocity of the 
flow closer to the sonic value than the low specific 
speed type. 

Since, in determining the flow velocity C,, use must be 
made of a contraction factor (defined as actual inlet area 
allowing for blade thickness divided by theoretical inlet 
area without allowance for any blade thickness), to 
account for the thickness of the vane inlet edges when 
calculating the impeller inlet area, it is important, if the 
impeller is a high specific speed type operating with a 
large flow, that the contraction factor be accurately 
determined and closely held to this figure in manufac- 
ture. 

If the inlet flow velocity C, is plotted against contrac- 
tion factor, as in Fig. 8, it will be seen that, if C, is 
initially high, slight manufacturing variations toward 
heaviness in inlet vane thickness may cause C, to rise to 
its sonic value and prevent the impeller from passing 
its designed flow. 


DERIVATION OF CURVES FOR AC 


In designing or testing a centrifugal compressor, the 
inlet conditions of pressure and temperature generally 
given are either those of the ambient air or of the air 
moving with low velocity in the inlet pipe and are, 
therefore, total pressure P, and total temperature T,. 
If the inlet flow velocity is calculated using these fig- 
ures, a value Cy will be obtained which does not take 
into account the change of density occurring as a re- 
sult of the acceleration of the air as it approaches and 
enters the impeller inlet. 

The true inlet velocity C, is given by C, = Cy) + AC, 
and AC may be quickly determined from curves de- 
veloped in the following manner:* From the well- 
known nozzle formula, 


Ci? = 2g J(Cp)(T, = Ti) = (1) 


where 


* The particular method of analysis used here was suggested by 
Mr. Jon Schenk. 


J = Mechanical Equivalent of Heat = 778 
ft.lbs. per B.t.u. 

Cp = Specific heat of air at constant pressure 

T; = Airstream, or static air temperature at the 
impeller inlet, °F. abs. 

AT,4. = Adiabatic temperature change from 7, to 


For an adiabatic expansion it can be shown that 
= T, (X/(X + 1] (2) 
where 


xX = | 

€p = the expansion ratio occurring during the ac- 
celeration of the air = P,/P, 

K = ratio of the specific heats = Cp/C, 


If an expansion efficiency of 95 per cent is assumed, 
based on the energy change, and the expression for A7,4, 
given in Eq. (2) is substituted in Eq. (1), then 


Ci? = 2gJCp(T,) [0.95X/(X + 1)] 
Putting 
2gJCp[0.95X/(X + 1)] = Ki 
and assuming that Cp is constant, then 


K, = 11495.0[X/(X + 1)] 


and 
CY = 
or 
T, = (3) 

Also, 

Ci = Of = Coler)(Ti/T,) 
and since 

AT = [0.95X/(X + LIT; 
then 
T, = T, {1 — [0.95X/(X + I]} 

or 


Ti/T, = 1 — [0.95X/(X + 1)] (5) 
and, substituting for 7;/T, from Eq. (5) in Eq. (4), 
= Coep{1 — [0.95X/(X + 1)]} 


Putting 
Kz = 1/ep{1 — [0.95X/(X + 1)]} 


Ci = Co X (1/Re) (6) 


then 


Assuming various expansion ratios ep, the correspond- 
ing values of X may be calculated, or looked up in 
existing tables, and a tabulation made along the lines 
shown in Table 1. 

Knowing ep, 1/K,; and 1/K2 and assuming different 
values of Cy, the corresponding figures for AC, Ci, and 
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TABLE 1 
X 0.95 X 1 1 0.95 X 5 
€; X X¥+41 X¥+41 +41 Ky X¥+1 
Pes 0. 02734 1.02734 0. 0266 0.0253 306 0.00327 0.9747 0.932 1.072 
1:2 0.05295 1.05295 0. 0503 0.0478 578 0.00173 0.9522 0.875 1.142 
1.3 0.07708 1.07708 0.0716 0.068 824 0.00121 0.932 0.825 1.292 
1.4 0.0999 1.0999 0.091 0. 0864 1,047 0.000955 0.9136 0.781 1.279 
TABLE 2 
ep = 1.30 C, = G X 1/K2 = 606 631 655 
1/K, = 0.00121 C2 = 367,000 398,000 429,000 
T, = Cj? X 1/K; = 444 481 519 
= 1.212 AC=G-GQG= 106 111 115 


7, may be obtained most conveniently in a tabulated 
form (see Table 2). 

Having tabulated the desired range of ep and Cy, the 
results may be plotted in the form of a set of curves as in 
Fig. 9. 

The uncorrected approach velocity Cy may be cal- 
culated from the given total temperature and total 
pressure conditions and the curves then used for ob- 
taining the correction value AC corresponding to that 
particular Cy value at the given total inlet tempera- 
ture. Then the corrected approach velocity C; = Cy + 
AC. 

The proper theoretical inlet vane angle may then be 
determined, as illustrated in Fig. 1, and, having de- 
cided on the angle of attack to be used, the actual vane 
inlet angle is known. 
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On Airfoil Theory and Experiment 


ALAN POPE* 
Georgia Institute of Technology 


SUMMARY 


The more exact equation for the slope of the lift curve is re- 
stated, and a relation for the location of the aerodynamic center 
is advanced. An approximate method for determining the mini- 
mum drag coefficient and extent of the low drag bucket is also 


given. 
LirT AND MOMENT 


wc THE EARLY YEARS Of airfoil theory, lift 
curve slopes were so far below the thin airfoil value 
of 2x per rad. that it became customary to use that 
criteria for thick airfoils as well, and both in writing 
and conversation the 27 value (0.1095 per deg.) became 
a standard that was broadly accepted. In a similar 
way the thin airfoil approximation, that the aerody- 
namic center was at the quarter-chord, became ‘“‘stand- 
ard,’ and, indeed, for a long period all airfoil tests re- 
ported demonstrated aerodynamic centers a small 
distance ahead of the quarter-chord. 

When it became possible in the light of newer in- 
formation to design airfoils with more advantageous 
pressure distributions, disturbing results appeared. 
Lift curve slopes were higher than the first approxima- 
tion theory, and aerodynamic centers appeared behind 
the airfoil quarter-chord. It became at once necessary 
to re-examine the approximations previously made for 
thin airfoils in order to include the effect of thickness. 

To begin with, let us note the equation for the lift 
curve slope as shown in reference 1 (with one addi- 
tional term of the series retained). 


[1 + (0.77d/t)] 
da [1 + (0.774/1)*] 


where d/t is the thickness ratio. ’ 

It is seen that the lift curve slope is a function of air- 
foil thickness, increasing with thickness, and should 
only be the famous 27 when the thickness is zero. 

We find the theoretical location of the aerodynamic 
center for a symmetrical airfoil as follows: 

From page 86 of reference 2, the moment of a simple 
airfoil about the center of its parent circle is 


M. = (p/2)V*4mc* sin 2a 
where c equals the distance from the origin to the circle’s 


rear stagnation point. Hence, the moment coefficient 
about the aerodynamic center will be 


Cmac = Cme — €(0.5 — a.c.) 
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where ¢, = lift coefficient and a.c. = chordwise location 
of the aerodynamic center. Hence, 
sin 2a 
[4e(1 + &)}? 
where = 0.77(d/t) and: 4c(1 + = airfoil chord. 
Next, 


Cmac = 2(1 + — ¢(0.5 — ac.) 


c(0.5 — a.c.) 


and applying the condition that démac/dc,; = 0, we get, 
finally, 
ac. = [1/411 + &)(1 + €)] — '/2 
Thus, the theoretical location of the aerodynamic 
center is also a function of airfoil thickness, being at the 


26.7 per cent point of an airfoil 9 per cent thick and at 
the 27.6 per cent point for one 15 per cent thick. 
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Fic. 3. Demonstration of how the effective airfoil slope changes 
with angle of attack. 


Comparisons of theory and experiment for both lift 
curve slope and aerodynamic center position are shown 
in Figs. 1 and 2. The experimental points are from 
reference 3. 

It is well known, and quite apparent in the figures, 
that the perfect fluid values of lift curve slope and aero- 
dynamic center position are not attained in practice. 
The answer is that the above theories have no pro- 
vision for viscosity. Further, its effect is not symmet- 
rical about the airfoil mean line so that the so-called 
effective airfoil, about which the nonviscous theoretical 
flow would realize, is changed as the airfoil angle of 
attack and, hence, the boundary layer are changed. 
Pinkerton‘ made a striking explanation of the situation 
when he reduced the Theodorsen pressure distributions 
so that they agreed more closely with those actually 
experienced with the early type airfoils. The steps 
necessary for this procedure could be interpreted as 
the same steps that one would follow if the camber of 
the airfoil were decreasing: with angle of attack. Re- 
stated, the boundary layer on the upper surface, pro- 
ceeding against an adverse gradient more severe than 
that on the lower surface, becomes thicker than that of 
the lower surface and effectively adds S-camber with 
each angle of attack. Thus one could say that the lift 
curve we plot is a composite of a great number of theo- 
retical airfoils, each developing the theoretical lift 
curve slope but having decreasing camber and, hence 
being moved to the right on the figure. 

For example, consider a 65-012 airfoil. At a = 0, 
we have c,; = 0. At a = 1.0, we should have, accord- 
ing to the theory, c, = 0.119. However, at 1.0°, the in- 
creased boundary-layer thickness on the upper surface 
has produced an effective S-camber that changes 
az... from 0° to +0.075, and, hence, the value of c, for 
1° becomes 0.119(1.0 — 0.075) = 0.110. This is demon- 
strated in Fig. 3. 

We see that the effect on the lift curve slope is tied 
in with the difference in rate of change of boundary- 
layer thickness on upper and lower surfaces and would 


reach the theoretical slope only if upper and lower 
boundary layers grew at the same rate. 

Indeed, since to develop lift we must have a lower 
pressure and, hence, a greater adverse pressure gradient 
toward the trailing edge on the upper surface, it is hard 
to see how the boundary layer could ever add more 
thickness to the lower surface than to the upper and 
ever yield a slope greater than the theory. As far as 
the author has been able to ascertain, no slopes higher 
than the theory have occurred. 

One frequent variation is noted: A change in regu- 
larity in the pressure distributions that ia turn varies 
the rate of boundary layer build-up so that breaks in 
the lift curve occur and several slopes develop. 

The effect on the moment coefficient is of equal in- 
terest. The boundary-layer S-camber should decrease 
the moment coefficient with angle of attack. The 
usual approximation that the moment is due entirely 
to the lift enables us to make a simple presentation of 
this fact. If the values of ¢n (1/4) are plotted against 
¢,, the location of the aerodynamic center may easily 
be shown to be 


a.c. = 0.25 — 


or, if the moment curve has a positive slope, the aero- 
dynamic center is ahead of the quarter-chord. If it 
has a negative slope, the aerodynamic center is behind 
the quarter-chord. Since the aerodynamic center 
position and slope of the lift curve are functions of the 
S-camber, we would expect an airfoil with low lift 
curve slope to have a forward a.c. as great S-camber is 
indicated. On the other hand, airfoils whose pressure 
distributions encourage extended thin boundary layers 
have less camber change and more nearly approach the 
theoretical lift curve slope and aerodynamic center 
position. The highest lift curve slope is noted for the 
63000 series, which also has the most rearward a.c. 
This is strong corroboration for Pinkerton’s explanation. 

In closing this section on lift and aerodynamic center, 
it is in order to note that, while limiting values have 
been indicated for lift curve slope and aerodynamic 
center position, no attempt has been made to predict 
how closely a particular new airfoil will approach the 
limit. In general, it can only be said that airfoils whose 
minimum pressure peaks occur well back from the 
leading edge will have steeper lift curve slopes and 
more rearward aerodynamic center positions than the 
older forward pressure peak airfoils. 


DRAG 


Because of the pronounced effect of several variables, 
a prediction of the complete drag curve of a new airfoil 
is probably not possible. But, fortunately, in the low 
drag region, which is of considerable interest, a short 
approximate method can be employed that for modern 
airfoils of 9 to 18 per cent thickness ratio yields the 
drag ‘‘bucket’’ shape and extent and the value of the 
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drag coefficient within usual engineering limits of ac- 
curacy. First of all, it is in order to review the sources 
of profile drag and the behavior of the boundary layer 
in order to have a complete picture of the problem of 
drag prediction. Usually the profile drag is divided into 
skin friction and pressure drag. In the region under 
discussion (the low-drag range), the method to be pre- 
sented takes the pressure drag as zero, and any error in- 
curred thereby is included in the prediction of the skin 
friction. As a matter of interest it has been noted 
elsewhere,® and been corroborated by the author, that 
the true pressure drag is most difficult to obtain. Plots 
of the chordwise pressure distribution at low angles of 
attack consisting as they do of the small difference in 
large areas are not easily evaluated. 

The skin friction is due directly to the viscous forces 
in the boundary layer. These may be of three types: 
laminar from front stagnation point to a little aft of the 
minimum pressure point; transitory for the next few 
per cent chord; and, finally, turbulent over the re- 
mainder. Formulas have been advanced for the drag of 
laminar and turbulent boundary layers in regions of 
constant pressure, but the studies of the drag of transi- 
tory flow and on the effect of the pressure gradients on 
the constant pressufe drag values have not been too 
successful. For an approximate drag value, however 
the author has found that good agreement exists be- 
tween the drag calculated using the skin friction 
fo mulas for constant pressure regions and that actually 
measured in the wind tunnel under fairly low turbu- 
lence conditions. 

The procedure is as follows: 

(1) Plot the theoretical pressure distributions of the 
airfoil in question according to the method of reference 
3or6. Every degree from zero lift to a = 8 is adequate. 

(2) Examine the pressure patterns. 

(3) Assume that laminar flow will exist up to the 
minimum pressure points shown on upper and lower 
surfaces. For the desired speed ‘and airfoil size, calcu- 
late the Reynolds Number of the laminar flow and its 
drag coefficient from 


Cai = 1.328, ‘s/RN 


Assume that turbulent boundary layer will exist over 
the remainder of the surface. However, it should have 
somewhat less drag than over a flat plate of similar 
chord, since it starts with the thickness of the laminar 
boundary layer. The assumption is made that the drag 
coefficient should be based on the Reynolds Number of 
the turbuent region plus one-half the Reynolds Number 
of the laminar region, using the formula 


Cat = 0.455/ (logy RN)?* 
The drag coefficients must then be reduced to apply 
to their proper wiag area, of course. 


This whole process may best be examined by an 
example: 


Fic. 4. 


Suppose the pressure distribution of a 2-ft. chord air- 
foil indicates laminar flow to exist to the 40 per cent 
chord point on both upper and lower surfaces. Find the 
drag coefficient at 200 ft. per sec under standard condi- 
tions (p/u = 6,380); turbulence free conditions. 

(1) The Reynolds Number of the laminar flow will be 


RN, = (p/n) V1 = 6380 X 200 X (0.4 X 2) = 1,020,000 
and for both surfaces will be 
Ca, = 2(1.328/+/1,020,000) = 0.00263 


(2) The Reynolds number of the turbulent region 
will be based on a chord of 0,6 X 2 + 1/,(0.4 X 2) or 
1.6 ft. 


For both surfaces 
0.455 
(logio 2,080,000)" 00772 


(3) Now to adjust the coefficients so that they apply 
to the complete airfoil cz, we must take 0.4c,, and 
0.6ca:. Hence, 


Ca = 0.4(0.00263) + 0.6(0.00772) = 0.0057 


The method outlined above has been applied to a 
modified 65-311 shown in Fig. 4. The comparison be- 
tween calculated and measured drag is shown in Fig. 5. 
The agreement appears to be within experimental 
error. 
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Fic. 5. Comparison of approximation and experimental drag 
values for an NACA 64-311 (modified) at RN. = 2,200,000. 
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As the complete comparison entails considerable 
labor, values of mininium drag only were examined for 
the 64, 65, and 66 series symmetrical airfoils. In every 
case, the skin friction calculated as above came within 
Ac, = 0.0008, and usually less than Ac, = 0.0004, of 
the value measured in the tunnel. These are presented 
in Figs. 5, 6, and 7. 


CONCLUSIONS 


The more exact formula for the slope of the lift curve 
has been restated, and a relation for the location of the 
aerodynamic center advanced. A short method for de- 
termining the approximate drag coefficient and extent 
of the laminar flow ‘‘bucket’’ is also presented. This 
drag approximation is most applicable to smooth low 
drag airfoils from 9 to 18 per cent thick at Reynolds 
Numbers from 1 to 7 million. 
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Dynamic Lateral Stability as Influenced by 
Mass Distribution 


LEONARD STERNFIELD* ann MARION O. McKINNEY, JR.* 


National Advisory Committee for Aeronautics 


SUMMARY 


The present paper presents a résumé of some recent theoretical 
and experimental investigations of the effects of mass distribution, 
wing incidence, wing loading, and altitude on the stability of the 
lateral oscillation of airplanes. The results are presented in the 
form of charts that show calculated oscillatory stability bound- 
aries and time histories of the lateral oscillations of models in 
flight in the Langley free-flight tunnel. 

On the basis of these lateral stability studies, the preparation 
of general stability charts for the prediction of the stability of the 
lateral oscillation does not appear feasible. Small variations in 
some of the mass and aerodynamic parameters may cause such 
a pronounced change in the oscillatory stability that a special 
stability analysis should be made for each airplane. 


INTRODUCTION 


— IN THE DESIGN of recent high-speed air- 
planes has been toward high wing loadings with 
either thin low aspect ratio wings or swept wings. 
The combination of these two factors make it necessary 
to distribute most of the weight of the airplane along 
the fuselage or to carry a large part of the fuel load in 
external tanks, which are often located at the wing 
tips. These design practices have led to unusually high 
rolling or yawing moments of inertia. Another design 
trend has been toward the use of relatively large angles 
of wing incidence to reduce the fuselage drag for high 
altitude or cruising flight or to reduce the fuselage 
ground angle and thereby simplify the landing-gear de- 
sign. In order to obtain the angle of attack of the wing 
desired for a given flight condition, the principal longi- 
tudinal axis of an airplane with positive wing incidence 
would be located below that of the airplane with no 
wing incidence. 

Numerous investigations to determine the effects of 
mass distribution on dynamic-lateral stability have 
been made at the Langley laboratory of the N.A.C.A. 
and more work is now in progress. The present paper 
is a résumé of the theoretical and experimental investi- 
gations of various phases of the problem that have 
already been completed. The effects of varying the 
wing loading, the rolling and yawing moments of 
inertia, and the inclination of the principal longitudinal 
axis of inertia relative to the flight-path axis are dis- 
cussed. 


Presented at the Aircraft Design Session, Sixteenth Annual 
Meeting, I.A.S., New York, January 26-29, 1948. 

* Stabiltty Research Division, Langley Memorial Aeronautical 
Laboratory. 


The results of the various mass-distribution studies 
are not presented in detail. The principal findings are 
illustrated by specific examples and reference is made to 
several N.A.C.A. publications on the subject for more 
general presentations. The results are presented herein 
in the form of calculated oscillatory-stability boundaries 
or time histories of the lateral oscillations of models 
flying in the Langley free-flight tunnel. 


SYMBOLS AND COEFFICIENTS 


The forces and moments acting on the airplane are referred to 
the stability axes, which are defined as an orthogonal system of 
axes intersecting at the airplane center of gravity in which the Z 
axis is in the plane of symmetry and perpendicular to the relative 
wind, the XY axis is in the plane of symmetry and perpendicular 
to the Z axis, and the J) axis is perpendicular to the plane of sym- 
metfy. A diagram of these axes showing the positive direction of 
forces and moments is presented in Fig. 1. 

The symbols and coeff.cients are defined as follows: 


o = angle of bank, rad. except where otherwise noted 

p = rolling angular velocity, rad. per sec. (d@/dt) 

vy . = angle of azimuth, rad. except where ctherwise noted 
Y 


Wind 
direction 


direction 


Fic. 1. The stability system of axes. Arrows indicate positive 
directions of moments, forces, and control surface deflections. 
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= yawing angular velocity, rad. per sec. (dy /dt) 

= angle of sideslip, rad. except where otherwise noted 
(v/V) 

= angle of incidence of wing, deg. 

= sideslip velocity along the Y axis 

= airspeed, ft. per sec. 

= mass density of air, slugs per cu. ft. 

= dynamic pressure, lbs. per sq. ft. (!/29 V?) 

= wing span, ft. 

= wing area, sq.ft. 

= weight of airplane, lbs. 

= mass of airplane, slugs (W/g) 

= acceleration of gravity, ft. per sec. per sec. 

=" relative-density factor (m/pSb) 

= angle of attack of principal longitudinal axis of air- 
plane, positive when forward end of principal axis 
is above the X axis 

= angle of flight path to horizontal axis, positive in a 
climb, deg. 

= radius of gyration in roll about principal longitudinal 
axis, ft. 

= nondimensional radius of gyration in roll about prin- 
cipal longitudinal axis (kx,/b) 

= radius of gyration in yaw about principal normal 
axis, ft. a 

= nondimensional radius of gyration in yaw about prin- 
cipal normal axis (k z,/b) 

= nondimensional radius of gyration in roll about longi- 
tudinal stability axis cos? 7 + sift? 

= nondimensional radius of gyration in yaw about ver- 
tical stability axis ~/ Kz? cos? 7 + Kz? sin? n 

= nondimensional product-of-inertia parameter 
[—(Kzo? — sin cos 7] 

= trim lift coefficient (W cos y/qS) 

= rolling moment coefficient (rolling moment/gSb) 

= yawing moment coefficient (yawing moment /gSb) 

= lateral force coefficient (lateral force/qS) 

= effective dihedral derivative, rate of change of rolling- 

moment coefficient with angle of sideslip, per rad. in 

equations and per deg. in figures (0C;/08) 


= directional-stability derivative, rate of change of 
yawing-moment coefficient with angle of sideslip, 
per rad. in equations and per deg. in figures (O0C,/ 
op) 

= lateral-force derivative, rate of change of lateral-force 
coefficient with angle of sideslip, per rad. (OCy/08) 

= damping-in-yaw derivative, rate of change of yawing- 
moment coefficient with yawing-angular-velocity 
factor, per rad. [0C,/0(rb/2V) ] 

= rate of change of yawing-moment coefficient with 
rolling-angular-velocity factor, per rad. [0C,/0(pb/ 
2V)) 

= damping-in-roll derivative, rate of change of rolling- 
moment coefficient with rolling-angular-velocity 
factor, per rad. [0C,/0(pb/2V) | 

= rate of change of rolling-moment coefficient with 
-yawing-angular-velocity factor, per rad. [0C;/ 
O(rb/2V)] 

= rate of change of lateral-force coefficient with rolling- 
angular-velocity factor, per rad. [0Cy/0(pb/2V) | 

= rate of change of iateral-force coefficient with yawing- 
angular-velocity factor, per rad. [OCy/0(rb/2V)] 

= time, sec. 

= distance in spans ( Vt/b) 

= differential operator (d/ds) 


EQUATIONS OF MOTION 


The method of analysis of dynamic stability is based 
on the classical theory of Bryan,' who was the first to 
apply the theory of small oscillations to the dynamics 
of mechanical flight. 

The nondimensional linearized equations of motion, 
referred to the stability axes, used to calculate the 
oscillatory-stability boundary for any flight condition 
are: 


Roll: 
(2u(Kx*D*@ — = Ci,8 + *, 2€, Do + 
1/2C,Dy 
Yaw: 
2u(Kz?D*y — KxzD*o) = C18 + + 
Sideslip: 


2u(DB + Dy) = Cyg8 + '/2Cy,Do + Cio + 
VeCy,Dy + (Cz tan 


It can be shown? that, when ge’ is substituted for 
we for and for in the equations written in 
determinant form, \ must be a root of the equation 


+ + AEX + F=0 
where 


A= Kxz”) 

B= + Kx?Cy, + 
2K x2*Cyg + KxzCi, + KxzCry) 

C = + 4uKx*Cag + Kz*CipCr, + 
+ KxzCi,Cyg + 4uKyzCig + 
K7?Cy,Cig Kx?Cy,Cng KxzCy,Cig) 

k= 7 + + 
7 2uCiKxzCug 2uC Kz?Cig 
Crg C, tan y — yzCigCy tan y + 
Cyp + 

F = — + '/2C, tan y X 
(CipCrg CrpCig) 


The conditions necessary for neutral oscillatory sta- 
bility are that the coefficients A, B, C, and E be positive 
and Routh’s discriminant, R = BCE — AE? — B?F, 
be equal zero. The oscillatory stability boundary is 
plotted as a function of the directional-stability deriva- 
tive C,,, and the effective-dihedral derivative C,z. For 
an airplane designed with a combination of C,, and 
Cig located above the boundary, oscillatory stability 
will occur. 


RESULTS AND DISCUSSION 


Effect of Product of Inertia 


Until recently, the product-of-inertia effect has usu- 
ally been neglected in lateral-stability analyses because 
the lateral-stability study of reference 3 indicated that 
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this effect was relatively unimportant for conventional 
airplanes. It has been shown in references 4, 5, and 6, 
however, that the product of inertia may have a pro- 
nounced effect on the lateral-stability of present-day 
airplanes designed for high-speed flight because of high 
wing loadings, large differences between rolling and 
yawing moments of inertia, high operational altitudes 
and the aerodynamic characteristics of low aspect ratio 
or swept wings. 

The product of inertia results from the inclination of 
the principal longitudinal axis of inertia relative to the 
reference axis which, in lateral-stability equations 
based on the stability axes, is the flight path. This 
angle of attack of the principal axis causes the inertia 
forces to produce a coupling between the rolling and 
yawing motions so that a rolling acceleration produces 
a yawing moment and a yawing acceleration produces 
arolling moment. The product of inertia factor 


Kxz = —(Kz,? — Kx,?) sin 1 cos 9 


appears in the rolling and yawing moment equations 
and is a function of the mass distribution and the angle 
of attack of the principal longitudinal axis of inertia. 

In order to illustrate the importance of the product 
of inertia on lateral stability, calculated oscillatory sta- 
bility boundaries are presented in Fig. 2 for a high- 
speed airplane with a wing loading of 70 lbs. per sq-ft., 
cruising at an altitude of 30,000 ft. The boundaries are 
plotted as a function of the directional-stability deriva- 
tives C,, and the effective-dihedral derivative Cig for 
two cases that represent the inclination of the principal 
axis at an angle of 2° below the flight path at the nose 
(n = —2°) and no inclination of the principal axis to 
the flight path (7 = 0°). A comparison of the two 
boundaries shows a large destabilizing shift in the 
boundary as the principal axis falls below the flight 
path. That is, the airplane is oscillatorily stable only 
for those combinations of C,, and Cig located above 
the boundary, and, therefore, as the boundary shifts 
upward, the stability at any particular combination of 
C,g and Cy, will be reduced. 

The sideslipping motion of the airplane subsequent 
to a disturbance in sideslip is shown in Fig. 3 for two 
values of 7, —2° and 0°. The airplane was assumed 
to have values of C,, and Cyg of 0.0045 and —0.0025, 
respectively (point indicated on Fig. 2). Although 
the motion is stable for both cases, the oscillation is 
less damped when the forward end of the principal axis 
is inclined downward with respect to the flight path. 
When the principal axis is aligned with the flight path, 
the oscillation damps to one-half amplitude in 2.50 
sec.; whereas, when the principal axis is inclined below 
the flight path, the oscillation requires 4.75 sec. to 
damp to one-half amplitude. The period of the oscilla- 
tion is about 1.5 sec. for both cases. Hence, the num- 
ber of cycles required for the oscillation to damp to half 
amplitude is 1.67 and 3.17 for principal-axis inclinations 
7 of 0° and —2°, respectively. 
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Effect of the angle of attack o: the principal longitudinal 
axis on the oscillatory-stability boundary. 
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The damping of the oscillation would increase still 
further if the angle of attack of the principal axis n were 
increased in a positive direction. Calculations of sta- 
bility boundaries for various values of 9 have indi- 
cated, however, that the most pronounced effect of 
product of inertia on the stability boundary occurs at 
small values of ». This fact is clearly illustrated in 
Fig. 4, which shows a cross plot taken from conven- 
tional stability boundaries in the C,,Ci, plane to pro- 
duce a stability boundary as a function of 9 and 
C,g for a constant value of Cig. For any particular 
value of n, stability will exist for values of C,, above this 
boundary. It is apparent from this chart that a given 
incremental change in 7 will produce a more pronounced 
change in the stable range of values of C,, when 7 is 
small than when 7 is large. 

Several high-speed airplanes have recently been de- 
signed with relatively large angles of incidence of the 
wing relative to the fuselage in order to reduce the 
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Calculated time history of the motion in sideslip for a 
hypothetical high-speed airplane. 


Fic. 3. 


.005 
.003 
.002 
1 for 
| 
; 
ive 
” 
1S 5 
va- 
or ‘ 
: 
4 
. 
4] 
2 3 5 6 
u- 
se 


AERONAUTICAL SCIENCES-—JULY, 


414 JOURNAL OF THE 
DIRECTIONAL 
STABILITY 
DERIVATIVE, 
n 
008 - 
Cy) =-0.002 
0064 
.004- 
.002- 
-.002- 
2 6 8 10 
7, DEG 
Fic. 4. Effect of the angle of attack of the principal longitudinal 


axis on the minimum value of C,g required for oscillatory stability. 
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fuselage drag for high-speed flight at high altitude or 
reduce the ground angle and thereby simplify the land- 
ing-gear design. With positive wing incidence, however, 
the airplane will be less stable than with no incidence. 
This result is caused by the fact that, for a given angle 
of attack of the wing, the principal axis for an airplane 
with positive wing incidence is located below that of the 
airplane with no wing incidence. 


In order to illustrate this effect of product of inertia 
on lateral stability, a flying model was tested in the 
Langley free-flight tunnel with the wing set at two 
angles of incidence, 0° and 10°.® In each of these 
configurations the model was flown at the same lift 
coefficient that corresponded to an angle of attack of 
10° for the wing. Thus, the angle of attack of the 
principal axis that coincided with the fuselage reference 
axis was 10° when the wing was set at 0° incidence and 
was 0° when the wing was set at 10° incidence. The 
calculated stability boundaries for the model in these 
two configurations are shown in Fig. 5 along with a 
test point that locates the model on the chart according 
to its values of C,,, and Cig. The boundaries indicated 
stability for the model when the incidence was 0° and 
instability: when the incidence was 10°. This fact was 
verified by flight tests of the model. Fig. 6 shows meas- 
ured time histories of the uncontrolled motions of the 
model for the two test configurations. For the sake 
of simplicity, only the rolling motions are presented. 
The period and damping of the yawing and sideslipping 
motions, however, were identical with those of the roll- 
ing motion. The time histories show a damped 6scilla- 
tion when the wing incidence was 0° and an oscillation 
of increasing amplitude when the wing incidence was 
10°. Thus, these data show a qualitative check on the 
calculated effects of wing incidence and product of 
inertia. The more complete data presented in reference 


. 6 also show a good quantitative check of the stability 


calculations. 


The general flight behavior of the models was found 
to be influenced to a pronounced degree by the stability 
of the lateral oscillation. The two model configurations 
were identical except for the wing incidence, but gen- 
eral flight behavior was satisfactory when the lateral 
oscillation was stable (i,, = 0°) and definitely unsatis- 
factory when the lateral oscillation was unstable 
(i,, = 10°). 


Effect of Mass Distributed Along the Fuselage 


The effect of varying the yawing moment of inertia 
by changing the mass distribution along the fuselage is 
illustrated in Fig. 7. This figure shows calculated oscil- 
latory stability boundaries for two loading conditions 
which correspond approximately to a conventional air- 
plane (kz, = 4.8) and a modern high-speed design (kz, = 
9.6) for two angles of attack of the principal axis, 7 = 
0° and n = 5°. 
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Variation of the yawing moment of inertia has little 
effect on the oscillatory stability boundary when the 
principal axis is at 5° positive angle of attack. This 
result can be explained by analyzing the effect of varying 
the yawing moment of inertia on the product of inertia. 
An increase in the yawing moment of inertia causes an 
increase in the magnitude of the product of inertia 


factor 


Kxz = —(Kz,? — Kx,”) sin 7 cos 9 


which has a stabilizing effect on the oscillatory stability. 
As the moment of inertia in yaw is increased, however, 
the airplane finds it more difficult to yaw, and the reac- 
tion moment due to the product of inertia, caused by 
yawing acceleration, become less stabilizing. Thus, the 
two effects tend to oppose each other, and the net re- 
sult is an extremely small shift in the oscillatory sta- 
bility boundary due to variation of the yawing moment 
of inertia. - 

When the principal axis was aligned with the flight 
path (y = 0°) the results of the computations shown on 
Fig. 7 show an entirely different trend. As the yawing 
moment of inertia increases the oscillatory boundary 
shifts upward and causes a decrease in the stable region. 

A comparison of the results for the two cases (n = 0° 
and » = 5°) indicates that if the mass is distributed 
largely along the fuselage it is almost necessary that the 
principal axis be inclined above the flight path in all 
flight conditions in order to obtain oscillatory lateral 


stability. 


Effect of Mass Distributed Along the Wing 


In general, an increase in the rolling moment of in- 
ertia has been found to decrease the stability of the 
lateral oscillations of an airplane. There may be ex- 
ceptions to this general statement for certain combina- 
tions of values of the moment of inertia in roll and the 
angle of attack of the principal axis. The effects of 
varying the rolling moment of inertia are discussed in 
the theoretical analysis of reference 5 and the experi- 
mental study of reference 7. 

For illustration in the present paper, Fig. 8 shows the 
effect of varying the rolling of inertia on the oscillatory 
stability boundary of a model tested in the free-flight 
tunnel. The rolling moment of inertia was increased 
by moving a given weight from the wing center section 
out to the wing tips to simulate a change from a belly 
fuel tank to tip tanks. This means of varying the roll- 
ing moment of inertia obviously involves a simultane- 
ous change in the yawing moment of inertia. This case 
is considered a practical one, however, inasmuch as 
changes in the rolling moment of inertia generally in- 
volve variation of the loading along the wing and con- 
sequent changes in the yawing moment of inertia. 

As shown in Fig. 8, an increase in the rolling moment 
of inertia causes a large shift in the oscillatory stability 
boundary in such a direction as to indicate a decrease 
in stability. The boundaries indicate that for the com- 
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Fic. 7. Effect of radius of gyration in yaw on the oscillatory- 
stability boundary. 


bination of C,,3; and C)g for which the model was tested, 
shown on Fig. 8 as a test point, the model was stable 
when the load was located at the center section and un- 
stable when the load was located at the wing tips. 
This change from stability to instability was verified 
by flight tests of the model as shown in Fig. 9, which 
presents time histories of the rolling motions of the 
model for the two conditions represented by the sta- 
bility boundaries of Fig. 8. These flight records clearly 
show stability for the low-moment-of-inertia condition 
and instability for the high-moment-of-inertia condi- 
tion. 

The general flight behavior of the model was judged 
to be satisfactory in the stable, low-moment-of-inertia 
condition. For the high-moment-of-inertia condition, 
however, the general flight behavior was definitely un- 
satisfactory. In fact, sustained flights were impossible 
because the model was so unstable that the free-flight- 
tunnel pilot could not control it. This poor behavior is 
illustrated in Fig. 9, which shows the amplitude of the 
rolling motion to be increasing in spite of the pilot’s 
efforts to control it. 
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Fic. 8. Effect on the oscillatory stability of the change in 
mass distribution caused by moving the load from the wing center 


section to the wing tips. 
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Fic. 10. Effect of the relative density factor, » on the oscilla- 
tory-stability boundary. 


Effect of Wing Loading and Altitude 


An increase in the wing loading or altitude of an 
airplane generally has been found to be detrimental to 
its lateral stability. This result has been obtained from 
other investigations as well as from those given here in 
references 4 and 8. 

The effect of varying wing loading and altitude is 
illustrated in Fig. 10, which is a stability chart for a 
model tested in the Langley free-flight tunnel to study 
the effects of wing loading and altitude. The effects of 
these two parameters are treated simultaneously by 
considering variations in -the relative-density factor yu 
(ratio of airplane density to air density) because this 
factor varies directly with both wing loading and alti- 
tude (u = m/pSb). Fig. 10 shows oscillatory stability 
boundaries for various values of the relative density fac- 
tor. In order to provide orientation on this chart, 
Table 1 is presented. 

It is apparent from Fig. 10 that increasing wing load- 
ing or altitude causes a shift of the stability boundaries 


TABLE 1 


Wing Loading Altitude 


Airplane Type (Lbs. per Sq.Ft.)  (Ft.) 
Light plane 7 0 2.6 
10,000 3.5 
World War II fighter 40 0 13.0 
40 35,000 43.0 
Postwar high-speed design 100 0 60.0 
100 50,000 400.0 


in such a direction as to indicate a decrease in stability. 
These stability boundaries indicate that there is a pro- 
nounced effect of wing loading and altitude on stability 
for values of » less than 30, whereas these factors have 
little effect on stability when yu exceeds a value of about 
30. 

An experimental study of the effects of wing loading 
and altitude on the lateral stability of a flying model in 
the free-flight tunnel verified the pronounced effects of 
these factors for the range of values of » between 5 and 
30, which was the range covered in the tests. It was 
not possible to include higher values of » in these tests 
because the low air speed of the tunnel limited the 
maximum wing loading of the model. 

The general flight behavior of the model was found to 
be dependent to a considerable degree upon the oscil- 
latory stability. Thus, increasing » caused the general 
flight behavior of the model to become worse over the 
entire range of values of C,, and Cig covered in the 
tests. 


GENERAL REMARKS 


The present paper indicates in general the effect of 
mass distribution, wing incidence, wing loading, and 
altitude on the oscillatory stability. The results are 
illustrated for an airplane or model with a given set of 
values of mass and aerodynamic parameters. How- 
ever, as shown in more complete lateral stability stud- 
ies, small variations in some of the mass and aerody- 
namic parameters may cause a pronounced change in 
the oscillatory stability. On the basis of these detailed 
studies, therefore, the preparation of general stability 
charts for the prediction of the stability of the lateral 
oscillation does not appear feasible. It appears 
necessary to make a separate stability analysis for each 
airplane. 


CONCLUSIONS 


From the results of the present paper, it may be con- 
cluded that: 

(1) Small angles of inclination of the principal longi- 
tudinal axis of inertia relative to the flight path cause 
a pronounced increase in the oscillatory stability. 
Further increases in the angle of attack of the principal 
axis, however, do not result in a proportionate increase 
in the oscillatory stability. 
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DYNAMIC LATERAL STABILITY 


(2) An increase in wing incidence results in a reduc- 
tion in oscillatory lateral stability. 

(3) An increase in the rolling or yawing moments of 
inertia, either independently or simultaneously, gener- 
ally appears to cause a reduction.in oscillatory lateral 
stability. 

(4) An increase in wing loading or altitude generally 
causes an appreciable reduction in lateral stability. 
For the ease of high wing loadings and altitudes, how- 
ever, there is virtually no effect of varying these fac- 
tors on the oscillatory stability. 
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The Rolling Moment Due to Sideslip 
for a Swept Wing 


JOHN W. MILES* 
Unwersity of California at Los Angeles 


SUMMARY 


Using strip theory, the rolling moment coefficient due to side- 
slip (Cig) is computed for a swept wing, a result that has been er- 
roneously derived in the literature.'!:? Multhopp’s theory* for a 
wing in sideslip and a modification of Prandtl’s elliptic wing 
theory are used to correct the result for induction effects. The 
effects of geometrical dihedral, twist, camber, elevons, or flaps, 
and fuselage* are also treated. The results are compared with 
available experimental data,°~’ and it is found that, if Weiss- 
inger’s* value of Cig for a straight wing is added to the result fora 
swépt wing, the agreement is satisfactory for most purposes. 
It is remarked that these results are subject to modification by 
aeroelastic deflection of the wing. 


NOTATION 

b = wing span, tip to tip 

C1 = section lift coefficient 

Cm, = section pitching moment coefficient at c; = 0 

c = wing chord, measured parallel to air stream 

d = outboard distance of center of pressure from root 
chord, expressed as a fraction of the semispan 

mo = section lift curve slope in the absence of sweep 

v = sideslip velocity, positive to starboard 

q = dynamic pressure = (1/2)pU? 

= longitudinal axis, positive in direction of undisturbed 

flight 

y = lateral axis, positive to starboard 

A.R. = aspect ratio = b?/S 

CL. = lift coefficient = L/gS 

C; = rolling moment coefficient = (rolling moment) /(g.Sd) 

Cy = pitching moment coefficient = (pitching moment) 
(qSc,,), where c,, is the mean aerodynamic chord, 
and the moment is evaluated about aerodynamic 
center of mean aerodynamic chord 

Cu, = pitching moment coefficient at C; = 0 

L = lift force : 

ly = lift force of unswept wing 

A = total wing surface area 

U = free-stream velocity 

a = angle of attack measured as a rotation about (y) axis 

perpendicular to air-frame plane of symmetry 

a = angle of attack of section at c; = 0 ; 

B = sideslip angle = (v/L’) 

6 = elevon deflection, measured as a rotation about hinge 
line, positive for down elevon 

€ = structural twist, measured as a 


y = angle of yaw, positive for starboard wing displaced 
aft 


Received January 16, 1948. 
* Assistant Professor of Engineering, Department of Engi- 
neering. 


p = mass density of air 

o = sweepback angle (negative for sweep forward) of*line 
of aerodynamic centers 

T = ratio of tip chord to root chord 

Tr = dihedral angle, measured as a rotation about wind 


axis, positive for wing-up; this angle is infinitesimal, 
so that it is immaterial whether [ is measured as 
a rotation about the wind axis or about a body axis 
that is at an angle of attack a with respect to the 
wind axis 


INTRODUCTION 


— EFFECT OF SWEEP on the rolling moment due 
to sideslip has been treated by Betz,! Weissinger,? 
and many others in the literature. As Weissinger points 
out,” almost all of the treatments are incorrect, and, 
while he avoids an error of his predecessors, he him- 
self makes a different error. While the following 
treatment is only approximate in many places, it 
attempts to clarify some of the errors in previous 
analyses. 


Betz’ TREATMENT 


Betz! argues (incorrectly) that the lift on a swept 
wing is reduced by cos? o, where o is the sweepback 
angle. Hence, if 8 is the sideslip angle (v/U) and Ly 
is the lift that the same wing panel would have if it 
were not swept, the actual lift of the swept wing may be 
written 


L = Ly) cos*o (1) 


and the incremental lift (positive for the starboard wing 
when v is to starboard) for the starboard panel alone 
is given by 


0 
= — 2 = i 
2 cos? (¢ — B) sin cosa 


2 dB 
Ltano (2) 


Then, if d is the spanwise location of the center of pres- 
sure expressed as a fraction of the semispan, the roll- 
ing moment coefficient per unit sideslip angle is given 


by 


b 
~ O&(v/U) q Sb 


Cig = —dC, tan o (3) 


where the full span (5) is utilized in defining the 
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moment coefficient, and the reference axis is the direc- 
tion of undisturbed flight. * 


Lirt COEFFICIENT FOR A SWEPT WING 


Actually, as Weissinger points out,’ the lift on a 
swept wing is reduced by cos a instead of cos* o, pro- 
vided that the unswept and swept wings being com- 
pared have the same angte of attack measured as a 
rotation about an axis perpendicular to the air-frame 
plane of symmetry. It is the looseness in stating the 
axis of rotation which leads to the confusion in speci- 
fying cos? o, for if the rotation is made about the swept 
axis (as in the case of an aileron), the lift coefficient is 
actually reduced by cos? c. 

Weissinger’s justification for the use of cos o is ob- 
scure, and the following argument is believed to be 
clearer. The free-stream velocity (U) may be repre- 
sented by its spanwise (U sin a) and normal (U cos a) 
components. If the angle of attack (a) is assumed 
small, it may be imagined to consist of a rotation (a@ cos 
«) about the swept axis, plus a component (—a sin o) 
about the axis parallel to (U sin a). The effective 
dynamic pressure is '/.p(U cos o)*, and it produces 
a component of lift proportional to a cos o and a second 
component proportional to the upwash produced by the 
spanwise flow and the dihedral; hence, the total lift 
per unit span is given by 


(—U sin (—asin (4) 
a — 
(U cos 2? da 


acosa (4) 


and the effect of sweepback on the section lift coefficient 
may be represented by writing 


(dc,/da)(c) = cos (dce,/da)(0) (5) 


It should be emphasized that a is measured as a rota- 
tion about the (y) axis perpendicular to the air-frame 
plane of symmetry, and in those cases where the rota- 
tion is about the swept axis, such as produced by aileron 
deflections, structural twist, cambered sections, etc., 
the correct reduction is cos’*¢. It may also be remarked 
that a rotation 6 about the swept axis introduces an 
effective dihedral component 6 sin ¢; the case of the 
elevon is treated as an example in the latter portion of 
the present paper. 


* In recent correspondence with the author, Dr. F. H. Clauser 
(The Johns Hopkins University) has suggested that Betz’ argu- 
ments presumed that the wing was displaced in such a fashion 
that the line of aerodynamic centers remains in the plane of the 
free-stream velocity at the angle of attack under consideration. 
In any event, his results are not directly applicable to the case 
at hand. 


CoRRECT DERIVATION OF Cig 


In view of the foregoing argument, Eq. (2) must be 
replaced by 


10) 
L 
tane (6) 
and Eq. (3) is replaced by 
Cis = (-—d '2\Cz tan (7) 


An alternative derivation of Eq. (7) is giveninthe Appen- 
dix. Weissinger’s result® for C,, (for a rectangular wing) 
differs from Eq. (7) by a factor of 2, caused by use of 
the semispan instead of the span in defining the rolling 
moment coefficient, and by stating sin o in place of tan 
ao. The latter discrepancy is an error on Weissinger’s 
part caused by confusing the lift coefficients of the 
swept and unswept wings. 


ASPECT RATIO CORRECTION 


The foregoing results are two-dimensional in the 
sense that they assume that the additional, asym- 
metric lift distribution resulting from sideslip pro- 
duces no changes in the vortex pattern. 

Actually, a wing having an asymmetric lift distribu- 
tion behaves, with respect to the effects of shed vor- 
tices, approximately as if the two panels were inde- 
pendent. (This approximation is exact in treating the 
rolling moment due to an asymmetric distribution over 
an elliptic plan form.) Moreover, the Prandtl theory 
for an unswept elliptic wing suggests that the two- 
dimensional slope of the lift curve for a swept wing 
should be reduced by aspect ratio approximately in 
accordance with the relation 


dC, My COS 
—(AR.,o) = — (8) 
da 1 + (my cos o/n A-R.) 

my = (dC,/da) (~, 0) (9) 


Weissinger apparently uses this line of argument in 

suggesting that the reduction of the lift coefficient for 

a finite wing due to sweepback is given approximately 
A.R. + 2 ) dC, 


by 

dC, ( 

— (AR., = | — (ALR., = 
a 4) A.R.seca+2/ da 


A.R. dC, 
+ ;) 4 de (A.R.,0) (10) 


which relation checks satisfactorily with more exact 
calculations for wings having sweepbacks up to 45°.” 

Now, in accordance with Eq. (8), if C, in Eq. (7) 
refers to the C, due to the wing attitude in the absence 
of the sideslipping motion being considered, the result 
for the additional lift due to sideslip should be cor- 
rected (approximately) by a factor 


(2) 
| 
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Mo COS o 2my cos 
where mp is the section slope of the lift curve in the 
absence of sweepback, as given by Eq. (9). 

On the other hand, Multhopp* has shown that the 
effect of the free vortex sheet (caused by the steady 
motion) on the sideslipping wing is to produce an addi- 
tional effective angle of attack distribution 


Aa = (11) 


where a; is the induced downwash resulting from steady 
flow, the positive and negative signs referring to the 
starboard and port wings, respectively, and 8 is positive 
for sideslip to starboard. Multiplying Aa by the slope 
of the lift curve, as given by Eq. (8), and assuming that 
the downwash is constant and equal to that which would 
be present for an elliptic plan form—namely, 


tano 


a = AR. (12) 
yields 
My COS o Cr 
AC, = { | t 1 
1 +(2m cos o/n A.R.)J AR? ms 


Calculating the rolling moment coefficient correspond- 
ing to Eq. (13) and adding the result to that of Eq. (7), 
as corrected for aspect ratio effects, yields 


* 11+ (m cos AR] 


d Mo COS CL 
2L1 + (2m cos A.R.) | ALR. 
d 
mea C,tane (14) 


so that the last two effects treated cancel one another. 
Weissinger? introduces the second correction, but not 
the first, so that his result is more complex. 


GEOMETRICAL DIHEDRAL 


If the geometrical dihedral of the wing is Ty, measured 
as a rotation about the wind axis,* then it is well known! 
that the rolling caused by sideslip is given by 


d (dC, 
Cup = 


If (dC,/da) is measured for symmetrical loading, it 
must be corrected for the smaller effective aspect ratio, 
in accordance with Eq. (8). Making this correction 
and adding the result to Eq. (14) yields 


(15) 


2 E 4 cos da 
A.R. 


(16) 


* See definition of ! under Notation. 
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In the case of a twisted or cambered wing, the di- 
hedral will exhibit a spanwise variation. Multhopp* 
assumes that the upwash due to dihedral should be 
measured at the three-quarter chord line and derives 
the following expression for the spanwise distribution 
of dihedral: 


dy dy 
2 dy dy dy 


where (x./4, %/1) are the coordinates of the quarter 
chord line (x measured positive forward in the direc- 
tion of undisturbed flight, y positive to starboard and 
perpendicular to the plane of symmetry, and z com- 
pleting a right-handed triad, positive down), €(y) 


represents the angle of structural twist (measured . 


from the zero lift angle for the entire wing), and 
a represents the local section of attack for zero (sec- 
tion) lift. 

Actually, the effective upwash acts at the three- 
quarter chord line only if the camber is parabolic (or 
lacking entirely),? and Multhopp’s result, as stated in 
Eq. (17), is therefore strictly correct only when ao(y) is 
caused by a parabolic camber. Nevertheless, it may 
be expected that the approximate use of Eq. (17) in 
the case of nonparabolic cambers will be sufficiently 
accurate for most practical calculations. In particular, 
the last three terms will generally be small, especially 
if the contribution of the elevons is treated as in Eq. 
(24). 

For a swept wing, by definition, 

= —tano (18) 


and the rolling moment and pitching moment due to 
the second term in Eq. (17)—i.e., due .to twist—are 
given by 


Cm 
Cig 


b 
2tano 
ao) c(y) ydy (19) 


which result may also be derived independent of Mul- 
thopp’s concept of dihedral. 

In most practical cases, then, it will be sufficiently 
accurate to add the result of Eq. (19) to that of Eq. (16), 
although in the case of wings having variable dz,/4/dy or 
high camber, it might be preferable to calculate the 
additional C,, directly from Eq. (17). 


CENTER OF PRESSURE LOCATION 


The center of pressure location d should be obtained 
for asymmetric loading. While such a calculation 
could be carried out by the Weissinger method,’ the 
labor is scarcely justified. If the center of pressure for 
symmetrical loading is found, either by calculation or 
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ROLLING MOMENT 


TABLE 1 
ALR. d—Weissinger d—Eq. (20) d—Eq. (21) 
0 5 1 0.439 0.5 0.424 
15 0.450 0.5 0.424 
30 0.463 0.5 0.424 
rr 0.481 0.5 0.424 
0 5 1/2 0.424 0.445 0.424 
13 0.433 0.445 0.424 
30 0.443 0.445 0.424 
rts 0.457 0.445 0.424 
0 6 1/9 0.424 0.445 0.424 
30 0.447 0.445 0.424 
0 10 1/2 0.427 0.445 0.424 
5 0.472 0.445 0.424 


by wind-tunnel tests, it should be sufficiently accurate 
in conjunction with Eq. (16). 

A simple calculation using strip theory shows that 
the center of pressure of a linearly tapered, untwisted, 
uncambered wing is given by 


d= ak (20) 
3 1+ (Ctip/Croot) 
while for an elliptic wing 
d = 4/3x (21) 


Weissinger’ has calculated the symmetrical lift distribu- 
tion for several swept wings and obtained the center 
of pressure locations (see Table 1), which are com- 
pared with the result predicted by Eqs. (20) and (21). 
While the error in the strip theory result of Eq. (20) 
is generally too large for the computation of pitching 
moments (insofar as trim and static stability are con- 
cerned), it is probably smaller than other uncertainties 
in the computation of C,g. For an unswept wing, it 
is evident that the assumption of an elliptic lift dis- 
tribution is more accurate than strip theory. 


EFFECT OF SYMMETRICAL ELEVON DEFLECTION 


In many cases (e.g., flying wings) sweepback is util- 
ized in order to give ailerons a lever arm in pitch and 
thus take the place of an elevator. The ailerons may 
then be termed ‘‘elevons.”’ If the elevon deflection 6 
is produced by a rotation about the swept axis, as is 
customary, the effectiveness is reduced by cos? o—viz., 


= (O0C,/06)(0) cos? 


In addition, the elevon deflection gives an effective 
component of dihedral 


Té = é6sine (23) 


Following the procedure used in obtaining Eq. (16), 


it is found that the contribution of a symmetrical elevon 
deflection is given by 


Cig = bs tan o 


where (OC,/06,) is the rolling moment coefficient per 
unit asymmetrical elevon deflection, 5s is the symmetri- 
cal elevon deflection required to, trim the airplane in 
neutral flight, and 5s and 6, are both positive for star- 
board elevon down. (Since 0C,/06, will be negative, 


(24) 


Cig will be negative for down elevons.) The value of 
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(0C,/06,) to be used in Eq. (24) is that caused by lift 
alone and does not include the contribution of cambe 
introduced by the elevons. 

Eq. (24) may, of course, also be applied to flaps, re- 
membering that (OC,/06,) is to be computed (or meas- 
ured) for asymmetrical flap deflection. 


FURTHER EFFECTS OF CAMBER 


The effect of camber on the dihedral distribution has 
already received attention. In addition to this effect, 
the section pitching moment coefficients due to camber 
will also affect the value of C),. 

If it is assumed that the camber is described about 
the swept axis, the moment (positive stalling) about 
the swept axis is given by 


dM = (+ isina +jcosa) X 
[(1/2)p(U cos cos sec = 


emo (+ tana + 7) (25) 


where C,,0 is the section moment coefficient due to the 
camber of a section perpendicular to the swept line of 
aerodynamic centers, (7, j) are unit vectors in the posi- 
tive (x, y) directions, and the top and bottom signs 
refer to the starboard and port wings, respectively. 

In steady flight, the pitching moment coefficient 
due to camber is then given by Eq. (25) as 


Curso = (2 cost ¢/CyS) c(y)emo(y)dy (26) 


while the contributions of the two wing halves to the 
rolling moment cancel. 

For a wing in sideslip the effective normal velocity 
component in Eq. (25) is U cos o = v sin o, instead of 
U cos o, and the contribution to the rolling moment 
coefficient is then given by 


0(dM)/0(u/U) = (—i sine + j cos a) X 
(pU? sin cos a)(€ COS sec (27) 
whence the rolling moment coefficient is found to be 


Cig —2(Cm/b) Crp tan® (23) 


EFFECT OF DRAG 


If the c.g. of the airplane is not on the wind axis, then 
the drag of the airfoil will also contribute to the rolling 
moment due to sideslip. If h is the elevation of the 
c.g. above the wind axis, the rolling moment about an 
axis through the c.g. and parallel to the wind axis, due 
to sideslip and drag, is given approximately by 


h 
Cis = 
My COS 
2my cos *) 
1+( A.R. 


where Cp is the total drag coefficient for the wing. 


oC, 
(c. tano + ag (29) 
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Eq., (29) is derived similarly to Eq. (16), assuming 
that the distribution of drag over the wing is similar 
to the distribution of lift. Since it is unlikely that the 
magnitude of Cig given by Eq. (29) will be of any 
practical importance, the latter assumption appears 
satisfactory. 


Cig FOR STRAIGHT WIN. 
In another paper,‘ Weissinger shows that a straight 
wing in sideslip develops a rolling moment given 
by 


x (0.15 + ] 
Ci Feal 0.06 (0) 


for a tapered wing, and 


x (31) 


—@5 


AR. - 
where x is an empirical constant for which he recom- 
mends the value 1.2. Inasmuch as this contribution is 
relatively independent of sweep, it should be added to 
those components of C)g due to sweep. 


OVERALL RESULTS 


Superimposing all the effects treated,.it is found 
that the net rolling derivative is given approximately 


by 
2m cos — 
A.R. ) 
1.2 (0.15 + 28, | 
— 0.05 2 
Feat 0.05 + 2tanoe X 


(6 6s + b - b Go tan’o (32) 


where C,,) is that part of C,,o due to twist, C,,* is that 
part due to camber, including the effects of both wing 
curvature and elevon or flap deflection, and C,, the 
lift coefficient, does not include the contribution of 
elevons. 


My COS 


d 
Cig tan o 


THE YAWED WING 


If a wing is placed at an angle of yaw (y, positive for 
starboard wing displaced aft), the effects are similar 
to those of sideslip (¥ = —), with the additional 
factor that the pitching moment vector is rotated with 
respect to the roll axis (which has been taken in the 
direction of undisturbed flight); hence, the two roll 
derivatives are related by 


Cy —Cig — (€m/b)Cu (33) 


Multhopp® includes this contribution of the pitching 
moment in C,s, but this appears to be incorrect, at 
least for the present wind axis. 
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EXPERIMENTAL MEASUREMENTS 


The derivative C,, is generally measured in the wind 
tunnel by yawing the model. 

Betz! cites some early measurements of Blenk® on a 
100 by 20 cm. rectangular wing with a G 387 profile and 
a sweepback of 30°. The results compared with those 
predicted by Eqs. (14) and (32) (assuming d = '/,) are: 


a Cy C,—Ea. (8) Cig—Ea. (14) Cig—Ea. (32) 
0.89 0.295 —0.21 —0.28 
7.8" 1.24 1.240 —0.89 —0.98 


so that the agreement is only in order of magnitude, 
It should be remarked that much of the error might be 
in the computation of C; and in the uncertainty as to the 
value of Cy, which has not been included in comparing 
Cy and Cu. 

More recently, Maggin and Bennett® of the N.A.C.A, 
conducted measurements on five different swept wings 
with the results shown in Table 2. Maggin and Ben- 
nett do not state whether their values of C,3 have been 
corrected for Cy,. 


TABLE 2 
ALR. Cr Cig Cig— Ea. (14) Cig—Ea. (32) 
42° 5.9 0.500 0.5 -0.17 —0.098 —0.112 
42° 3.0 0.707 0.5 -0.14 —0.116 —0.179 
42° 2.0 0.793 0.5 —0.21 —0.118 —0.231 
—38° 5.9 0.500 0.5 +0.13 +0.085 +0.071 
—38° 3.0 0.707 0.5 +0.11 +0.092 +0.029 


Theil and Weissinger’ conducted measurements on a 
straight wing and a similar swept wing with the follow- 
ing results: 


o A.R. T Cy Cw Cig—Ea. (14) Cig—Ea. (32) 
0° 5 1/2 0.5 +0.018 0 —0.021 
35° 1/3 0.5 +0.100 —0.078 —0 099 


While the agreement between Theil and Weissinger’s 
measurements and theory is excellent, it must be re- 
garded as somewhat coincidental in view of the discrep- 
ancies between theory and the measurements of 
Maggin and Bennett. On the other hand, the latter 
measurements do not seem to be entirely consistent, and 
the results for the wing of aspect ratio 3 are particularly 
anomalous. In any event, it appears that the agree- 
ment between theory and experiment is satisfactory, 
at least for the sweptback plan forms, since there are 
many other factors that will affect the overall value of 
Cig for a complete airplane. 

Theil and Weissinger also conducted measurements 
of Cy for the same wings with flaps deflected. They 
state, transposing Multhopp’s conclusion with an error 
in sign, that (2/A.R.)Cyo should be added to Cy, [the 
2/A.R., instead of (c,,/b), is caused by using the mean 
geometric chord for a lever arm reference in defining 
Cyo and the semispan in defining C,,], and their experi- 
mental results bear out their conclusion as to sign and 
order of magnitude. Actually, it appears that the 
effect of flap deflection is described by the third and 
fifth terms in Eq. (32), but since it is not clear from 
their report what values must be ascribed to (OC,/06) 
and Cy because of flap deflection, a quantitative check 
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js not possible. In any event, their use of (2/A.R.)Cyso 
jsalmost certainly incorrect. 


EFFECTS OF FUSELAGE AND TAIL 


The foregoing results apply only. to an isolated wing 
or horizontal stabilizer. In addition, there will be 
contributions from the vertical stabilizer and from inter- 
action between the fuselage and the wing and between 
the fuselage and the tail, although the present treat- 
ment is concerned only with the wing. 

Multhopp* has computed the additional C,, due to an 
elliptic fuselage mounted on an elliptic wing. For a 
swept wing of arbitrary plan form, his result may be 
modified to read 


_ 8) (te + b 
= + (2m) cos b 


hp + sin hy 


b/), 2 
= [(4/2) — 2n(2,/b)], Zp > hy/2 (34b) 
= [—(m/2) — 2r(zr/b)], < —hrp (34c) 


where h, and 6, are the heighth and breadth of the fuse- 
lage cross section measured at the 0.75 chord of the 
wing and zy is the distance of the root chord below the 
center of the fuselage. 

No experimental data is available to check Eq. (33) 
for a swept wing, but Multhopp cites certain measure- 
ments that verify the result for an unswept wing. 


AEFROELASTIC EFFECTS 


It should be remarked that the foregoing results are 
computed for a rigid wing and are subject to consider- 
able modification if the wing experiences elastic de- 
flection due to its own air loads. Assuming that the 
line of the aerodynamic centers is forward of the elastic 
axis, torsion tends to increase both angle of attack and 
dihedral angle, while. bending decreases the angle of 
attack and increases the dihedral angle. For a wing 
having appreciable sweep (say 30° or more), the effects 
of bending will predominate, with the result that the 
effective dihedral will be increased, although this will 
be partially compensated by an inboard shift of the 
wing loading, with a resultant decrease of the lever arm 
parameter d. In addition, aeroelastic twist introduces 
an effective Cyo. 


APPENDIX 


The following derivation of c,g is longer but more 
precise than that given in Eq. (6) and also demonstrates 
the proprietry of Eq. (15) for the contribution of geo- 
metrical dihedral. 


Fic. 1. Wind and rotation components. 


Consider first a swept wing at zero angle of attack 
lying in the xy plane. Now rotate the wing through 
an angle —T about the x axis, the angle I’ being defined 
as the dihedral. Next, bring the wing to an angle of 
attack a by rotating it through an angle a about the 
yaxis. If both aand TI are assumed to be infinitesimal, 
these two rotations are commutative and may be speci- 
fied by the vectors a and T indicated in Fig. 1. 

The total rotation can also be described by the com- 
ponents (with respect to x’, y’ axes) 


a 


acoso + TI sing 
(Al) 


I’ = -asine + T cosa 
as indicated in the figure. Now assume the airplane 
has a forward velocity along the x axis and is sideslip- 
ping along the y axis with a velocity BU when 8 is in- 
finitesimal. The velocity components along the x’, y’ 
axes are given by 


4 U(cos ¢ + 6 sin a) | 


v’ =U(—sino + 6 cos a) 


Il 


(A2) 


The velocity components of the wind with respect to 
the wing are indicated in the sketch. Now, if (dc,/da)o + 
is the section lift curve slope, uncorrected for sweep- 
back, the section lift is given by the lift caused by a’ plus 
the upwash caused by I'’v’—viz., 


dL/dy = (1/2)pu’*c(dc,?/da)y [a’ + (I''v’/u’)] 
Substituting Eqs. (Al) and (A2) in Eq. (A3) yields 
(strip theory) 
(dL/dy) = (1/2)pU°c(de,/da)y cos a X 

[(a + BY) (1 + B tan o)| 


(A3) 


(A4) 
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whence the section lift coefficient due to sideslip is given 
by 


Cig = = 0 = [(dc,/da)y cos a] X 
(a tang +T) = c¢, tana + (de,/da), T (AS) 


A third derivation has been suggested by Francis 
H. Clauser in a letter to the author. ‘I defined two 
successive small rotations [ and a. Next I defined a 
wind with a small angle of yaw (i.e., with components 
U cos 8, U sin 8, 0). Next, by analytic geometry, I 
found the direction cosines of the following: (a) the 
normal to the chord plane, (b) the line of aerodynamic 
centers, and (c) the wind. From these, it was easy, by 
means of the scalar product, to obtain the component 
of wind along the line of aerodynamic centers. This 
component of course gives no lift. Subtracting this 
vectorially from the wind itself gave the component 
of wind perpendicular to the line of aerodynamic cen- 
ters—i.e., the component effective in producing lift. 
The absolute magnitude of this component is used in ob- 
taining the effective dynamic pressure. The angle be- 
tween this component and the normal to the chord 
plane is the complement of the total angle of attack 
In carrying out this work, I found that on some of the 
components of the vectors it was necessary to carry 
higher than first order terms in I’, a, and 6 because of 
the canceling out later of the first order terms.”’ 
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A Note on Transverse Bending of Beams 
Having Both Translating and 
Rotating Mass Elements 


R. H. SCANLAN* 


Rensselaer Polytechnic Institute 


SUMMARY 


The usual equations of transverse vibration of beams take into 
consideration only the translation of the distributed mass of the 
beam and not its rotation. In particular, this is true for the 
equations applied to beams of variable mass and section (as, for 
example, aircraft wings), the structural characteristics of which 
are approximated by a matrix of deflection influence coefficients 
ata finite set of points. This note presents a simple extension 
of these ideas to include effects of inertia elements in the beam 
which rotate, as well as translate, during beam transverse vibra- 


tion. 


DISCUSSION 


M**’ MODERN AIRCRAFT WINGS have elements such 
as deep nacelles, tip tanks, turrets, or floats. 
During wing vibration these elements contribute bend- 
ing moments due to their rotation, as well as to their 
Often it is reasonable to neglect these 
To ascertain certainly 


translation. 
effects; occasionally it is not. 


the effect of rotatory elements, analysis should include 
them. 

The usual matrix representation for deflection y of a 
cantilever beam having (say four) arbitrarily selected 
sections to be used during analysis is 


6(1, 1) 6(1, 2) 6(1, 3) 6(1, 4) |] 
Veo? Ye} _ 6(2, 1) 6(2, 2) 6(2, 3) 6(2, 4) x 
Ys 5(3, 1) 6(3, 2) 6(3, 3) 4(3, 4) 
8(4, 1) 8(4, 2) 6(4, 3) 4) 
m,000] [1 
Om200 V2 
00m30 yz} (1) 
where y; is the transverse of deflection of beam section 


1, 6(7, 7) is the deflection of section 7 due to a unit trans- 
verse load at section j, m; is the mass assumed con- 
centrated at section 7, and w is the circular frequency of 
vibration. 


When rotatory effects are to be accounted for, this equation may be replaced by an analogous one that takes 
into account the additional inertia elements assumed to rotate about axes through and perpendicular to axis of the 


beam: 
br(1, 1) 2) 3) Sp(1, 4) 1) 2) 3) 4) | 
ye Or(2, 1) dp(2, 2) Sp(2, 3) Sp(2, 4) 1) 2) 3) 4) 
V3 1) bp(3, 2) 5p(3, 3) 4) 1) 2) 3) 4) 
L/w? br(4, 1) 5p(4, 2) 5p(4, 3) 4) 1) Sy (4, 2) 3) Sae(4, 4) x 
or(1, 1) or(1, 2) op(1, 3) or(1, 4) 1) oy(1, 2) 3) 4) 
ye" or(2, 1) op(2, 2) op(2, 3) op(2, 4) or(2, 1) 2) oy(2, 3) ou(2, 4) 
ys" op(3, 1) op(3, 2) op(3, 3) or(3, 4) ou (3, 1) 2) or¢(3, 3) or7(3, 4) 
ya’ | or(4, 1) op(4, 2) op(4, 3) ov(4, 4) 1) ou(4, 2) 3) ou(4, 4) 
Tm 0 0 0 0 0 0 OF Po | 
0 m0 0 0 0 0 0 Yo 
6.6 @ '@- M4 (2) 
000 0h 0 0 0 wh 
10 0 0 0 0 0 0 Ly’ 
where 5r(i, 7) = transverse deflection at section 7 due to a 
vi = transverse deflection of section 7 unit transverse force at section j 
yi = beam slope at section 7 
5(1, 7) = transverse deflection at section 7 due to a 
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unit bending moment at section j 
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or(z, 7) = beam slope at section 7 due to a unit trans- 
verse force at section j 


ou(i, 7) = beam slope at section 7 due to a unit bend- 
ing moment at section j 

mM, = mass assumed concentrated at section 7 

I; = polar mass moment of inertia of mass 
assumed concentrated at section 7 about 
a horizontal axis perpendicular to the 
beam elastic axis 

4,j = 1,2,3,4 

w = natural circular frequency of transverse 


beam vibration 


The values or, can be obtained by use of 
standard beam slope and deflection calculations for 
unit loads. In the usual elastic case [Eq. (1)], the 
influence coefficient matrix is symmetrical. This is 
also true of Eq. (2). This holds by virtue of the ex- 
tended form of the Reciprocal Theorem for elastic 
structures. Thus, if F; is a force applied at section 7 
and M; is a moment applied simultaneously at section 
ge 

FiMj6u(t, = 1) 
by the energy form of the Reciprocal Theorem; thus, 
as usual also 


J) = Sr(J, 1), oy(t, J) = 2) 


Hence, here too, the influence coefficient matrix js 
symmetrical. 

From this fact it follows that an orthogonality condi. 
tion can be obtained for the modal columns of Eqs. (2). 
This has the form 


4- + MsysV4 + + 
+ = 0 (4) 


where y, j’ represent deflection and slope in a mode 
different from that which yields y, y’. Thus Eqs. (2) 

‘ and (4) are completely analogous to the standard case 
concerned with beam deflections only. 

It appears that a structure divided into n segments 
for analysis by usual methods must now give rise toa 
(2n)* dynamic matrix rather than an n* one. However, 
it may be seen that a y’ equation need be carried for 
only those inertial elements having outstanding rota- 
tional contributions to the motion and not for all ele- 
ments. 

The usual extensions of Eq. (1) to free-free systems, 
coupled systems, etc., can be made immediately from 
Eqs. (2). Modes computed from Eqs. (2) contain the 
additional information of slope, as well as deflection, 
at the selected points. This serves to employ the addi- 
tional accuracy inherent in carrying extra matrix ele- 
ments. 

An example is presented to compare frequencies of a 
cantilever calculated first by assuming merely transla- 
tional elements and then later assuming that the tip 
mass has an inertial contribution that rotates the end 
of the beam. 


EXAMPLE 


Assume a weightless uniform cantilever beam of EJ = 103 X 10° Ibs. in.2 Let the beam be SO in. long with 
concentrated masses of 50, 75, 65, and 30 Ibs. assumed concentrated at points 20, 40, 60, and 80 in. from the canti- 


lever support, respectively. 


The influence coefficient matrix is then (for deflection in inches, rotation in radians) : 


1) Sp(1, 2) Sp(1, 3) Se(1, 4) 4) 2.59 6.48 10.35 14.23 0.194 
1) 6p(2, 2) 3) 4) 4) 6.48 20.70 36.20 51.80 0.776 

1) 6¢(3, 2) 3) 4) 4) = | 10.35 36.20 69.90 105.00 1.75 x 
5y(4, 1) 5p(4, 2) 3) 4) 4) 14.23 51.80 105.00 165.80 3.14 

opr(4, 1) op(4, 2) or(4, 3) 4) oy(4, 4) 0.194 0.776 1.75 3.14 0.0776 


The last mass of the beam (30 Ibs. at the tip) will be assumed to retain its spanwise position but will variously be 
assumed as: (1) located directly on the beam, (2) rigidly supported 10 in. below the beam, (3) rigidly supported 
30 in. below the beam. In the first case it contributes no rotatory inertia; in the others it has a rotatory moment of 
inertia of 3,000 and 27,000 Ibs. in.’, respectively. The equations of motion in matrix form are 


v1 2.59 6.48 10.45 

386.4 x 108} 2 6.48 20.70 36.20 
ys | = | 10.35 36.20 69.90 

vs 14.23 51.80 105.00 

ya! 0.194 0.776 1.75 


14.23 0.194 50 0 COO] 
51.80 0.776 07 0 00 2 
105.00 1.75 0 0 65 0 Of I ys 
165.80 3.14 0 0 O 30 Of I 
3.14 0.0776 0 0 0 ly 


where J = 0, 3,000 or 27,000 Ibs. in.* Solutions by iteration are given as follows: 


(Continued on page 434) 
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SUMMARY 
Since the problem of determining bending stresses due to tor- 
sion in a box beam is closely related to the problem of finding end 
moments in a continuous beam, it is possible to determine torsion- 
bending stresses by a procedure similar to the Cross method of 
moment distribution widely used in the analysis of continuous 
frames. 
The method, as applied to the simplest case—the prismatic 
box of rectangular section—is described in this paper. Cases in- 
volving flexible bulkheads and complete cutouts of one or both 
surfaces are considered. Expressions for constants analogous to 
the fixed end moments, stiffness, carry-over, and rigidity factors 
of moment distribution are presented, and application of the 
method is illustrated with numerical examples. 
NOTATION 
A = effective flange area, sq.in. 
a = bay length, in. 
b = web width, in. 
C = carry-over factor 
E = Young's modulus, lbs. per sq.in. 
F = fixed end flange load due to torque, lbs. 
f = fixed end flange load due to bulkhead distortion, Ibs. 
G = shear modulus, lbs. per sq.in. 
I = moment of inertia, in.* 
n = shear-flexure parameter 
a = axial load in flange, lbs. 
g = shear flow, lbs. per in. 
R = rigidity factor, lbs. per in. 
S = stiffness factor, lbs. 
r = torque, in. lbs. 
t = web effective thickness, in. 
é = deflection 
6,¢,¥ = rotation angles 
Subscripts 
b = bulkhead 
( = cover 
= secondary 
5 = spar 
r = due-to torque 
1,2,3,4 = numbers of flanges or webs 
INTRODUCTION 


Ws A BOX BEAM is loaded in torsion, secondary 
shear flows occur in the vicinity of concentrated 
torques and of structural discontinuities. These 
secondary, or torsion-bending, shear flows alter the 
customarily assumed 7°/2A» shear flow distribution in 
covers and spar webs and generate axial stresses in the 
beam flanges, often called torsion-bending stresses. 
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In wing boxes the resulting increments of shear flow 
and axial stress can be important, especially near the 
junction of open and closed bays and near fuselage 
intersections. 

The purpose of this paper is to describe a method of 
determining torsion-bending stresses in prismatic boxes 
of rectangular section—a cross section approximated by 
most wing boxes. 

Recent investigators of this problem have, in general, 
followed the pattern set forth in a comprehensive paper 
by Ebner.' Best known, in this country are the sim- 
plifications and extensions of Ebner’s work by Kuhn.? * 
Kuhn's assumptions are used in this paper to derive 
formulas for certain elastic constants to be used in 
applying the method. 

The method presented here is a solution by successive 
corrections without the use of simultaneous equa- 
tions. When bulkheads can be assumed rigid (as in 
most cases), the procedure for determining flange loads 
is similar to that used in determining end moments in 
a continuous beam using moment distribution. Cases 
involving bulkheads of relatively low rigidity require 
the introduction of flexibility corrections similar to the 
sidesway corrections used in applying moment distribu- 
tion to a frame having joints that move laterally. 


GENERAL CONSIDERATIONS 


To simplify the problem, the actual box cross sec- 
tions, with bending material distributed chordwise, are 
transformed as outlined by Kuhn‘ into the simplified 
cross sections shown in Fig. 1, in which the four corner 
flanges carry only axial loads, P, and the webs joining 
the flanges carry only shear flows, g. The area of each 
corner flange in this substitute or simplified box con- 
sists of the corner flange area in the actual box plus a 
certain fraction of the effective bending material dis- 
tributed between the corner flanges. The limiting 
maximum value of the fraction is '/5.° A reduced 
value should be used in the vicinity of structural or 
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Fic. 1. Simplified cross section of box. 
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Fic. 3. Equivalence of shear flow systems. 


loading discontinuities to account for shear lag effects. 
The effective web thicknesses of the simplified box are 
made equal to those of the actual box. 

The problem now becomes one of finding the flange 
axial loads, P, at bulkhead stations in the simplified box 
caused by torques applied at one or more bulkheads 
(Fig. 2). The box consists of four corner flanges of 
effective areas As, and A, (not necessarily 
equal) connected by four webs of effective thicknesses 
ty, ta, ts, and t, (which may also be unequal). Flange 
areas and web thicknesses are constant between adja- 
cent bulkheads. Bending rigidity of the flanges is as- 
sumed to be negligible. Bulkhead effective thicknesses, 
fy, can have any value other than zero. The following 
additional assumptions are made: 

(1) ‘The reaction opposing the loading is statically 
determinate, which means that external restraint 


against twisting can exist at not more than one bulk- 
head station. 
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(2) Axial movement of the flanges is prevented at 
the root bulkhead. 
(3) Internal shear flows in the spar webs and covers 
are constant between adjacent bulkhead stations. 
(4) Flange axial loads vary linearly between bulk. 
heads. 


SIGN CONVENTION AND NOTATION 


Positive forces, shear flows, and directions are as 
shown in Fig. 2. Signs of flange axial loads and web 
shear flows apply to those in flange No. 1 and web No. 1, 
respective'y. It should be noted that, using this con- 
vention, a positive shear flow in webs No. 1 and No. 2 
acts clockwise on the section inboard. 


STATICAL RELATIONSHIPS 


Any torque, 7, acting on the sides of the box as ver- 
tical and horizontal couples 7; and 7», respectively, can 
be considered as being equivalent to the two shear flow 
systems’ of Fig. 3: 

(1) A primary system, gr, statically equal to the 
torque, 7, and expressed by the familiar shell formula, 


gr = T/2Ao (1) 


(2) A secondary system statically equal to zero. 

The shear flow acting on any given bulkhead is the 
resultant of the external secondary shear flow (due to 
the external torque applied at the bulkhead) and that 
due to the internal shear flow (in the webs of the two 
adjacent bays). 

The internal secondary shear flow in either of the 
adjacent bays is a function of the unknown flange 
loads, P, and is therefore statically indeterminate. By 
considering a free body of one of the flanges, it is easily 
shown that its magnitude is given by 


where SP is the sum of the flange loads at either end of 
the bay. 

As further matters of statics, it is important to note 
that (1) the primary system, q7, prodices no flange 
loads and (2) the loads in the four flanges at any given 
cross section of the box are equal in (absolute) magni- 
tude and are equivalent to a doubly antisymmetrical 
group of four forces. 


DEFINITIONS AND SUMMARY OF FORMULAS 


Certain elastic constants must be computed before 
the distribution procedure can be applied. These are 
the “stiffness,” ‘‘carry-over factor,’’ “‘fixed end flange 
load,” and “rigidity.” These constants are defined and 
summarized below, first for closed boxes and then for 
open boxes. Since the means used to obtain the formu- 
las for these constants is not a part of the method of 
this paper, derivations are placed in the Appendix. 
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TORSION-BENDING STRESSES IN BOX BEAMS 


Closed Boxes 


Fig. 4 shows a free body of an isolated bay with axial 
movement of the flanges prevented at the far end, A, and 
with the box cross section at both ends prevented from 
distorting. On applying flange loads, P, the end cross 
sections of the spars undergo angular rotation of oppo- 
site sense. The force, Pg, at the near end required to 
produce a unit angle between the spars (9 = 1 in Fig .4) 
is the stiffness, S, of the bay. 


S = [Ebi/az(1/A)][(m + 4)/(n + 1)] (3) 


The ratio of the flange load at the far end to the 
flange load at the movable end is called the “‘carry-over 
factor,” C. 


C = —[(m — 2)/(n + 4)] (4) 


The fixed end flange load, F, is the flange load induced 
at either end of a bay because of an applied torque. 
It is assumed that the far end is built-in and that, at 
the twisted (near) end, axial movement of the flanges is 
prevented and the rectangular shape of the box cross 
section is maintained. 


F = qra[(m, — ms)/(n + 1)] (5) 


The rigidity, R, of a bay is the shear flow in the bay 
required to distort the cross section at the movable end 
through a unit angle (¢@ = 1 in Fig. 5). It is assumed 
that the shape of the box cross section is maintained at 
the far end and that axial movement of the flanges is 
prevented at both ends. 


R = —[B3Ebib2/a*Z(1/A)][1/(n + 1)] (6) 


It is convenient, though not essential to modify the 
formulas for S and R:to allow for the one end free con- 
dition that exists at the tip bay. When the stiffness 
is so modified, the carry-over factor at the inboard end 
of the tip bay equals zero. When the modified rigidity 
is used, flange loads at the outboard end of the tip bay 
caused by distorting the bay cross section are zero. 
Using basic principles of moment distribution it is 
easily shown that the modified stiffness, S”, at the near 
end of a bay with the far end free is given by 


S” = S(1 — C*) (7) 


and that the modified rigidity, R”, of a bay with one 
end free is given by 


R" = (R/2)(1 — C) (8) 
In the preceding formulas 


n = 


n, = 
n, = 


modulus of rigidity 


~ modulus of elasticity 
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Fic. 4. Stiffness of a bay. 


Rigidity of a bay. 


Fic. 5. 


equalsthe sum of the reciprocalsof the effective 
flange areas in the given bay cquals (1/A,) + (1/A2) + 
(1/A3) + (1/A4) 

= (be/te) + (bs/ts) and = + 
(b3/ts) equals the sum of the effective 6/t values of the 
covers and spars, respectively. 


Db/t = (b/t)e + (b/2), 


a = bay length 

by = bs = box depth 

bo = by = box width 

qr = primary shear flow in bay due to applied torque 
= T/2A9 


It is important to note that the flange load at either 
end induced by distorting the bay cross section through 
a unit angle equals aR. The flange load at the in- 
board end of a tip bay, regarded as a one end free bay, 
equals 2aR”. [See Eq. (2).] Signs of these flange 
loads are readily established by visualizing the elastic 
curve of the distorted bay. 


Open Boxes 

For boxes with one or both covers removed, definitions 
and notation similar to those for closed boxes apply. 
Expressions are as follows: 


S = Eb,/aX(1/A) (9) 
C= -1 (10) 
F = qra (11) 
R=0 (12) 
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Fic. 6. Constants for examples 1 and 2. 


For one end free: 


0 (13) 
R’ =0 (14) 


Bulkheads 


The rigidity, R,, of a bulkhead is defined as the shear 
flow required to deform the bulkhead through a unit 
shear angle and is given by 


R, = (15) 


where /, = effective thickness of the bulkhead. In 
the case of frame or truss-type bulkheads, t is the 
thickness of an “equivalent” solid web, which, under a 
given pair of shear couples, distorts through the same 
shear angle. 


APPLICATION 


The box shown in Fig. 6 loaded with a tip torque 
of 46.8 in. kips’ is used to illustrate the procedure. 
Computations have been made using a slide rule. 

Example | (Fig. 7) determines flange loads on the 
assumption that bulkheads do not distort. Usually, 
such distortions are unimportant, unless bulkhead 
rigidities are fairly small compared to the rigidities of 
adjacent bays. 

The effects of bulkhead flexibility are determined 
separately in example 2 (Figs. 8 and 9). 

Elastic constants are computed in Fig. 6 from Eqs. 
(3) through (15). 


The procedure for determining flange loads in Fig. 7 
is identical with that used in determining end moments 
in a continuous beam using the familiar routine of mo- 
ment distribution. The unbalanced flange load at 
each bulkhead station is distributed between adjacent 
bays in proportion to their relative stiffnesses, followed 
by carry-over of portions of the distributed flange loads 
to the far ends of adjacent bays. From the final flange 
loads obtained after six distributions, the secondary 
shear flows in the box are obtained using Eq. (2). 
Signs of the secondary shear flows are readily estab- 
lished by regarding the spar (or cover) as a beam loaded 
with end couples that are reacted by end shears. Total 
shear flows in web No. 1 (spar) and web No. 2 (cover) 
can be obtained as the algebraic sum of the primary 
(1'/2Ay) and secondary shear flows. Torsion-bending 
stresses can be obtained by dividing flange loads by 
flange areas. 

Figs. 8 and 9 determine the increments of flange load 
and shear flow in the box of Fig. 6 induced by bulkhead 
distortions, assuming that bulkhead A is infinitely 
rigid while those at B, C, and D have effective thick- 
nesses of 0.010 in. 

The auxiliary computations of Fig. 9 are used to 
maintain a running tabulation of the shear flows in the 
structure induced by bulkhead distortions and flange 
movements. The shear distribution factors (row 2) 
at any joint are proportional to the R values of the 
bulkhead and bays at that joint. At C and D, the 
modified value, Rep”, is used to allow for the actual 
free-end condition of bay CD. External secondary 
shear flows are recorded above the brackets in 
row 3. 

The solution of Fig. 7 assumes zero bulkhead distor- 
tion, which implies the existence of imaginary con- 
straints, supplying shear flows at bulkhead stations of 
sufficient magnitude to cancel the external and internal 
secondary shear flows that would otherwise distort the 
bulkhead. The net shear flow on the bulkhead is there- 
fore zero. These zero values and values of the internal 
secondary shear flows in the spars (from Fig. 7) are 

recorded in the first line below the brackets. 

On removing the imaginary constraint at any given 
bulkhead, the bulkhead and end cross sections of the 
adjacent box structure at the bulkhead distort through 
the same angle, each of these three elements absorbing 
a portion of the shear flow originally taken by the con- 
straint in proportion to the shear distribution factors. 
Thus at B, the total shear flow that comes on the 
structure because of removing the constraint is —52 — 
28.70 + 0) = —80.70. B therefore deflects downward 
and the —S80.70 is distributed between AB, bulkhead 
B, and BC in the proportions 0: 65.6: 34.4, giving the 
quantities 0, —52.90, and —27.80. A shear flow of 
27.80 is also induced at the far end, C, of bay BC, but 
for computing convenience these “far end’ values in 
this and any other distorted bay are not recorded. 
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Shear flows at C and D are distributed in similar 
fashion. 

The distributed shear flows just calculated induce 
fixed end flange loads, f, in the distorted bays, obtained 
as the algebraic difference between the distributed 
shear flow at the outboard end and that at the inboard 
end, multiplied by the bay length. (In the case of a tip 
bay regarded as a one end free bay, the multiplier is twice 
the bay length.) These algebraic differences are the 
values between the vertical arrows in Fig. 9. For ex- 
ample: 34.57 = 6.77 — (—27.80); 9.85 = 14.20 — 
4.35. The resulting fixed end flange loads are then cal- 
culated and recorded in Fig. 8. For example: 970 = 
34.57 XK 28; 551 = 9.85 (28 X 2), (one end free bay). 

Unbalanced flange loads at each bulkhead station are 
next distributed and carried over to the far ends of 
each bay. Flange loads are then determined by sum- 
ming the fixed end flange loads, distributed flange loads, 
and carry-overs, giving the values in row 7. (The re- 
sulting change in shear flows due to the change in flange 
loads from their original fixed end values is assumed to 
be absorbed by external constraints.) Shear flows in 
the box and bulkheads at this stage are as shown in 
Fig. 9, row 6. Shear flows to be distributed are again 
obtained as before. For example, at B: —52 — 
(—52.90 — 15.03 + 0) = 15.93 giving the distributed 
values 10.45 and 5.48. At C: — (15.03 + 12.95 + 
3.16) = —31.14 giving the values —8.75, — 16.75, and 
—5.64. At D: 52 — (—3.16 + 42.43) = 12.73 giving 
3.21 and 9.52. 

From these distributed shear flows, fixed end flange 
loads are again calculated and added to the previous 
summation (Fig. 8), giving a new set of unbalanced 
flange loads to be distributed and carried over. The 
procedure outlined is repeated to convergence, giving, 
after five cycles, the flange loads shown in Fig. 8. These, 
when added to those in Fig. 7, give the final values in- 
cluding the effects of bulkhead distortion. 


CONCLUDING REMARKS 


The procedure described herein enables a relatively 
simple determination of torsion-bending stresses, giving 
the same results as previously published methods with 
generally less effort. Its advantages increase with the 
number of bays considered. 

Two procedures are possible in cases involving flexible 
bulkheads. ‘That illustrated in examples 1 and 2 yields 
an initial set-of flange loads based on the assumption 
that bulkheads do not distort. A set of flange load 
increments due to bulkhead distortion only is then 
computed. Final loads are obtained by combining the 
two sets. One or two cycles will usually show whether 
the increments due to bulkhead flexibility are of suffi- 
cient importance to warrant the additional cycles 
needed to determine their true values. As an alternate 


procedure, a single-stage solution can be obtained by 
considering bulkhead distortions to occur from the 
start. Such a solution is identical in form to that of 
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example 2, with effects of bulkhead flexibility included 
in the resulting flange loads. 

More advanced torsion problems can be solved using 
the method of this paper. Among these are: (1) 
four flange prismatic boxes of trapezoidal cross section 
and (2) four flange tapered boxes of rectangular or 
trapezoidal cross section. To either of these con- 
figurations may be added combinations involving com- 
plete or partial surface cutouts and open bays having 
flanges with finite bending rigidity. 

Because of the additional variables in these cases, 
the formulas for the elastic constants become compli- 
cated. However, since the derivation of the formulas 
is not a part of the method of successive approximations, 
the distribution procedure itself is no more involved for 
these complex problems, provided that bulkheads can 
be assumed to be rigid. 


Appendix 


DERIVATIONS—CLOSED BOXES 


Fig. 2 represents a typical bay. Sections A and B are 
taken just inside of the bulkheads forming the ends of 


the bay. 

The forces acting are the shear flows qi, g2, gs, and q 
and flange loads P, and P,. It is assumed that section 
A is built-in. 

The conditions of statics are as follows: 

(A) Summation of forces = 0 giving q, = qs; and q = 
qs. (A positive shear flow acts clockwise when viewed 
from a point outboard of the given section.) 

(B) Internal torsion on section equals applied tor- 
sion: 


(q + = T 
from which 


= = — & (16) 


(C) Considering a free body of flange No. 1, assum- 
ing both P, and P, acting outboard (positive), gives 


Pa ga =P, + Ps qa = 0 
(17) 


From Eqs. (16) and (17) it is found that: 


qa = gra — 0.5P, — 0.5P, = qua (18) 
qa = gra + 0.5P4 + 0.5P2 = (19) 


Vertical deflection of spar No. 1 = 6, = bending de- 
flection plus shear deflection. 


= (a*b, (2P 4 + (qia/Gh) 


Substituting Eq. (19), noting that 1/J; = [(1/A1) + 
and rearranging, 


+ 


aG/ 1 1 


Sil 


La 


Simil 


= 


Spar 
equa 


pe 


in co! 


= 
Su 
1/A 
equa 
Stiffr 


luded 


using 
(1) 
ction 
at or 
con- 
com- 
aving 


‘ASES, 
mpli- 
1ulas 
ions, 
d for 
can 


TORSION-BENDING STRESSES 


Similarly, vertical deflection of spar No. 3 = 43: 


03 = t3 A A, 


a’G bs\ 
- — 2qra-¢ 


Lateral deflection of cover No. 2 = &:: 


b> = — (aby ‘6EI,)(2P 4 + (qoa ‘Gte) 


Substituting Eq. (18), noting that 1 J, = [(1 As) + 
| and rearranging, 


by a’G 1 
P,| —]) | — 2qra— 
Similarly, lateral deflection of cover No. 4 = 6,: 
1 { b,  2a°G 
bh @aG 
- — + —}|— 2g¢7a- 


Spar rotation = (4; — ¢,)/ = y,. Substituting the 
equations for 6; and 6;, noting that b:b. = dsbe, 


Cover rotation = (6; — 6,)/b; = y,. Substituting the 
equations for and noting that = d,b,, 


¢ equals change in shape of cross section equals change 
in corner angle = y, — y,. 


Stiffness and Carry-Over Factors 


o= 


From Eq. (20), setting ¢ = gr = 0, 


C = carry-over factor = P,/P,; = 
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1 
6, = > (sum of flange axial distortions) = 


by 
P,— l Ps, 
Substituting the value of P, from Eq. (4): 
t) + (a°G/: 3E) ~(1/A) 


t) + (4a2G BE) 3(1/A) 


S equals stiffness equals value of P, when 6, = 1. 


Eb, t Eby n+ 4 
A t A 
(3) 
Fixed End Flange Load 
From Eq. (20), setting ¢ = 0, P, = Px, and solving for 
Px, 


P, = fixed end flange load = F = 


2a 


Rigidity 
From Eq. (20) setting ¢ = 1, gr = 0, Pa = Pa = 
qa, 
q = rigidity = R = 
Gbyb be 1 
142 1 ) (6) 
n-+ 


DERIVATIONS FOR OPEN BOXES 
Assume that one or both covers of the box of Fig. 4 


are removed. 


Carry-Over and Stiffness Factors 


When flange loads P, are applied, statics requires that 
there be no shear flows in either spar. Therefore P, = 
Py». 


C = carry-over factor = P, ‘iP, = —1 (10) 


P, Ps, 


Stiffness = SS = value of Pg when 6 = 1 


S = Eb,/ad(1/A) (9) 


Fixed End Flange Load 


On applying torque, the spars deflect as beams with 
Because of symmetry, the inflection point 


fixed ends. 


4 

are 

Is of 

dq (5) 

tion 

wed 

16) 

1m- 

0 
17) 

8) 

9) 

b 2a?G 

> 

/ t 3E A 
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for axial load is midway between the ends From stat- 
ics, 
F = (2A0gr/be) (a/2)(1/bi) = gra (11) 
Rigidity 
Since the spars of the open box are assumed to have 
no resistance to deflections normal to their plane. 


R = rigidity = 0 (12) 


REFERENCES 


1 Ebner, H., Torsional Stresses in Box Beams with Cross Sec- 
tions Partially Restrained Against Warping, N.A.C.A. T.M. 
No. 744, May, 1934. 


2? Kuhn, P., Bending Stresses Due to Torsion in Cantilever Box 
Wings, N.A.C.A. Report No. 530, June, 1935. 

Kuhn, P., The Influence of Bulkhead Spacing on Bending 
Stresses Due to Torsion, N.A.C.A. Report No. 166, May, 1942. 

4 Kuhn, P., A Method of Calculating Bending Stresses Due to 
Torsion, N.A.C.A. Restricted Report No. 231, December, 1942. | 

5 Ibid., p. 28. 

6 Cox, H. L., On the Stressing of Polygonal Tubes with Particu- 
lar Reference to the Torsion of Tapered Tubes of Trapezoidal 
Section, British Air Ministry Reports and Memoranda No, 
1908, December, 1942. 

7 Kuhn, P., and Brilmyer, H. G., Stresses Near the Juncture of 
a Closed and an Open Torsion Box as Influenced by Bulkhead 
Flexibility, N.A.C.A. Wartime Report 91, originally issued 
August, 1945, as Advance Restricted Report L5G18. 


A Note on Transverse Bending of Beams 


Having Both Translating and Rotating Mass Elements 


(Continued from page 426) 


I=0 I = 3,000 Ibs. in.” I = 27,000 tbs. in.? 
Mode Mode Mode 
M 0.096 0.096 0.092 
y2 0.337 0.336 0.325 
¥3 0.655 0.654 0.642 
v1 1.000 1.000 1.000 
yy 0.0176 0.0184 
Frequency 570 cycles per min. 554 cycles per min. 536 cycles per min. 
Per cent change 0 3 6 
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Finite Deflections of Sandwich Plates 


ERIC REISSNER* 
Massachusetts Institute of Technology 


ABSTRACT 

This paper gives the basic differential equations for finite 
transverse deflections of sandwich plates under the following 
assumptions. The plate consists of a core layer and of two face 
layers of such construction that the face-parallel stresses in the 
core and the variation of the face stresses over the thickness of 
the face layers are negligible. The resultant equations permit 
the analysis of the effect of transverse shear stress deformation 
and transverse normal stress deformation in the core on the over- 
all behavior of the plate. It is shown that, in general, the effect 
of the transverse normal stresses in the core is negligibly small 
compared with the effect of the transverse shear stresses. The 
equations, when simplified by the omission of the transverse nor- 
mal stress terms, are brought into a form suitable for the solution 
of rectangular-plate problems. 

It is further shown that the range of deflections for which 
the linear (‘‘small deflection’’) theory is adequate decreases in 
accordance with a simple explicit formula as the core is made 
softer relative to the faces. 

The finite deflection equations are used to obtain plate-buck- 
ling equations that include the effect of transverse shear stress 


deformation on buckling loads. <A typical buckling problem is 


solved and discussed. 


(1) INTRODUCTION 


- THIS NOTE we derive a system of differential equa- 
tions for the small finite deflections of sandwich 
plates. This system of equations is a generalization of 
the now well-known results for homogeneous plates, 
where the problem can be reduced to two simultaneous 
equations for the transverse deflection w of the middle 
surface and for an Airy stress function F.! 

We consider a sandwich plate consisting of a core 
layer of thickness h — t and two face layers of thickness ¢ 
each. We assume that ¢ is small compared with h and 
that the values of the elastic constants E,, G; for the 
face layers are large compared with the values of the 
elastic constants F, G, for the core layer. We further 
assume that the products tE,, tG;, are large compared 
with the values of hE,, and hG,. 

On the basis of the assumption that t < h, we may 
assume that the stresses in the faces parallel to their 
planes are distributed uniformly over the thickness of 
the face layers.| On the basis of the assumption that 
hE, < tEy, we may neglect the face-parallel stresses in 
the core layer and their effect on the deformation of the 
composite plate. Thus, we treat the sandwich plate as 
a combination of two plates without bending stiffness 
(the face layers), and of a third plate (the core layer) 
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+ By this assumption we exclude the possibility of investigat- 
ing local (skin) instability. 


offering resistance only to transverse shear stresses and 
transverse normal stresses. 

A question of primary interest for such a composite 
plate is whether nonlinear effects have to be taken into 
account as soon as the transverse deflections are of the 
order of the face layer thickness ¢ or whether this is only 
the case when the transverse deflections become of the 
order of the core thickness 4. We shall decide this 
question on the basis of the system of equations ob- 
tained in this note. ' 

We somewhat limit the general applicability of the 
present results by the assumptions made concerning 
order of magnitude of thickness ratio and ratio of elastic 
constants. These assumptions should, however, be 
adequate in those cases where metal face layers are 
separated by a balsawood or synthetic material core 
layer. More general results can be obtained along the 
lines of the following work with the help of somewhat 
lengthier algebra. 

We restrict attention to face layer materials that 
are isotropic in the plane of the layers. This restric- 
tion also may be removed should the eventual usefulness 
of the present considerations warrant further develop- 
ments. 


(2) EguiLisrRIuM EQUATIONS AND STRESS-STRAIN 
RELATIONS FOR FACE AND CORE LAYERS 


The face layers are treated as thin plates with negli- 
gible bending stiffness normal to their planes. In the 
equilibrium equations for these plates (membranes), 
the effect of small finite displacements is taken into 
account. The loads applied to the membranes consist 
in external loads p and in the stresses at the junction of 
face layers and core layer. We inditate quantities 
referring to the upper membrane ( = +//2) by a sub- 
script + and quantities referring to the lower membrane 
(g = —h/2) by asubscript —. 

With the notation of Fig. 1 the equilibrium differ- 
ential equations for the membranes are the following: 


(ON, £/Ox) + (OS4/Oy) * tre = 0 (1) 
(OS ./dx) + (ON, = 0 (2) 
Ow 
2 4+2(s. + 
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Fic. 1. Element of sandwich plate. 


The equilibrium equations for the core layer are, 
under the assumption of negligible face-parallel core 
stresses, 


= 0,  dr,/dv = 0 (4, 5) 
(07,/dx) + (d7,/dy) + = 0 (6) 


The stress-strain relations for the face membranes 
are known to be of the following form: 


Ou « 1 OW < 1 

Ov 1 Ow 1 

Ou Ov Ow OW 1 


The stress-strain relations for the core layer may 
be written as 


= 0w/0z = a./E, (10) 
= (Ou/0z) + (Ow/Ox) = 7,/G, (11) 
Yy = (00/02) + (Ow/dy) = 7,/G, (12) 


The differential Eqs. (4) to (6) and Egs. (10) to (12) 
for the core layer must be integrated, and the results of 
the integration must be combined with the remaining 
equations for the face layers in such a way that a 
system of differential equations for the composite plate 
is obtained. 


(3) Tue DirrFERENTIAL EQUATIONS FOR THE Com- 
POSITE PLATE 


We begin with the integration of Eqs. (4) to (6) and 
Eqs. (10) to (12) for the core layer. From Eqs. (4) and 
(5) it follows that 7, and 7, do not vary across the thick- 
ness of the core. We define transverse shear stress 
resultants V, and V, by means of the following equa- 
tions 


V, = hr, (13, 14) 


V, = hr, 


1948 


From Eq. (6) it follows then, by integration, 
(OV,/0x) + (OV,/Oy) + = 0 (15) 


From Eq. (10) and from the fact that ¢, varies linearly 
over the thickness, it then follows that 

w+ — w- = + o,_)/2E, (16) 

Eqs. (11) and (12) imply the following relations: 

= [0( w dz) /Ox] + — u- 


Vy = [0( w dz) oy] + v+ — v- 


(17) 
(18) 
We may further write 
w dz = (Ow/dz)z dz = 
(w+ + w-) — (02/E,)2 dz 
As o-,is a linear function of z, we have 
= + + (22/h)(o.4 — 
and therewith 


Si wdz = 


h/2 


(19) 


(1/o)h(w+ + w-) — 


(h 2/ 1 2E,) z+ (20) 


We may take the last term in Eq. (20) from Eq. (15) and 


have then 
h? (OV, , OV, | 


fe) 
u+ — 


= aE (w, + w_) + 
h? (2V, , Wy | 


G, Ox 
v+ —v- (22) 


(21) 


h 
Ee 


Eqs. (16), (21), and (22) are the stress-strain re- 
lations for the core layer in a form suitable for use 
in the derivation of the equations for the composite 


plate. 
For the equations of the composite plate we define 
appropriate variables as follows: 


a= — u-)/h, (vy — 2-)/h (23) 


represent effective changes of slope of the normal to the 
undeformed middle surface; 


w = (we + w-) 
represents the effective transverse deflection of the 
middle surface; 

u = ('/2) (u+ + u-), v = (1/2) (vt + v-) 


represent the effective tangential components of dis- 
placements of the middle surface; and 


(24) 


(25) 


e = (w+ — w-)/h (26) 


represents the effective transverse normal strain for the 
composite plate. 


Finally 


for th 
layer, 
means 
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define 
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Fror 


Eqs. 
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_* 
? 
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(19) 


(20) 


») and 


(21) 


FINITE DEFLECTIONS 


In addition to the transverse shear stress resultants 
y, and V’, of Eqs. (13) and (14), we define stress result- 
ants and couples for the composite plate as follows: 


S=S+S- f 

(Nat —Ne-)/h, My = (Ny 
H = (S+ — S-)/h 


(27) 


Finally, we write 


o2— ) (29) 


= ("/2) + 


for the effective transverse normal stress in the core 
layer, and we define external load terms p and q by 
means of the following relations: 


P= = ("/2) (p+ — p-) 
Note, that once we have determined the quantities 
defined in Eqs. (23) to (30), we can revert to the original 
quantities referring to face layers and core layers indi- 
vidually through the use of Eqs. (23) to (30). 
We now obtain the equations of equilibrium for the 
composite plate by combining the six equations, Eqs. 
(1) to (3), by means of suitable addition and subtrac- 


tion. 
_From Eq. (1) follows: 


(30) 


(ON,/Ox) + {0S,/dy) = 0 (31) 
(0M,/Ox) + (OH /dy) — V, = (32) 
From Eq. (2) follows: 
(0S/dx) + (ON,/Oy) = 0 (33) 
(0H /dx) + (OM,/oy) — V, = 0 (34) 


From Eq. (3) are derived, after some transformations, 
the following two relations: 


Ny, 3 (a1 + 2H +M, 
=) =0 (35) 


Eqs. (31) and (33) are the usual equations of hori- 
zontal force equilibrium for elements of the plate. 
Eqs. (32) and (34) are the usual equations of moment 
equilibrium. 


Eq. (35) is the condition of transverse force equilib- 
tium. This equation contains terms that do not occur 
when homogeneous isotropic plates are considered— 
namely, the nonlinear terms that have the couples 
M,, M,, H and the transverse stress resultants V, and 
V, as factors. 
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Eq. (36) has no counterpart in the theory of homo- 
geneous plates, in the sense that the corresponding 
equation for the homogeneous plate contains informa- 
tion that is not of practical interest and is therefore 
never formulated. The significance of Eq. (36) is that 
it gives the local change of thickness of the plate caused 
directly by the external loads through the term g and 
resulting indirectly from the external loads by way of 
the nonlinear terms having stress resultants and couples 
as factors. 

We next obtain stress-strain relations by combining 
Eqs. (7) to (9) according to Eqs. (27) and (28). The 
resultant equations become, after some transforma- 
tions, 


dx \dx 4 UE, 


(37) 


4 N, — » Nz (38) 
Oy 2 Ox 
ou Ov wow h®0ede S (39) 
dy Ox | 4 2G, 
Oa Owde M, — v,M, 
ox dx Ox (1/2)th*E, 
Ow Oe My (41) 
ow re) H 
Oy Ox (1/2)th?G, 


In addition to the above six stress-strain relations, we 
have Eqs. (21) and (22) and Eq. (16), which may be 
written in the following form: 


(43 
hG, 12E, ox oy 

e= @,/E. (45) 


It is seen that in the relations involving N,, N,, S, 
M,, M,, and H there occur terms that are neglected for 
homogeneous isotropic plates—namely, the terms with 
e which represent the effect of transverse core flexibility. 
Eqs. (43) and (44) take into account the transverse 
shear deformability in the core. They reduce to what 
is known as Navier’s hypothesis when G, = E, = @. 

We shall show that the nonlinear terms involving e in 
Eq. (35) may be expected to be unimportant in all but 
the most extreme cases so far as softness of the core ma- 
terial is concerned and niay therefore, in general, be 
disregarded. 

We shall also show that the order of magnitude of the 
transverse deflections for which the relation between 
load and deflection becomes sensibly nonlinear depends 
on the relative softness of core and face material in a 
well-defined way. 
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(4) DIMENSIONAL ANALYSIS OF THE DIFFERENTIAL 
EQUATIONS FOR THE COMPOSITE PLATE 


To estimate the relative order of magnitude of the 
terms occurring in Eqs. (31) to (45) we may proceed 


as follows: We write 
N,, Ny, S & toa (46) 
M,, M,, H = tho, (47) 
Vi, Vy hr - (48) 
0/0x, 0/Ox = 1/1 (49) 


and introduce Eqs. (46) to (49) into equilibrium equa- 
tions and stress-strain relations. 

We write the resultant equations, with dimensionless 
functions c,, which are of order unity, in the following 
form: 


(tho,/l) + ahr ~ 0 (50) 
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+ + + (cathore/l?) ~0 (I 


Ee — + (cshtoae/I?) + (cetopw/l?) ~ 0 (59) 
)+ c3(w/l)? + co(he/l)? ~ oa/E, 


(C10 (Cu we /]?) =~ op hE, (54 
a+ (Cyw/l) & + (55 


We take 7 from Eq. (50) and og from Eq. (53) ang 
substitute in Eq. (51) in the form: 


ht 
p+ (1 + se) + 


/w\? w (eh\? 


We combine Eqs. (50), (54), and (55) and obtain 


ht Ey ht | 
=r 


[1 + PE, 


hw 
[2 


E, C19 (57) 


Combination of Eqs. (57) and (56) gives a relation of the form 


th® (w\? w en (h/1)? (ht 


A similar transformation changes Eq. (52) into 
| ht Ey 

ri 

l l 


1 a3 
+ Cos E, 

We restrict further attention to cases for which Eq. 
(59) amounts essentially to the relation 


q ht 


“p E. 


e~@q/E, (59a) 


This should be true in most cases of practical interest. 
We further assume that g ~ p and see that then the 
terms involving e in Eq. (58) are negligibly small, pro- 
vided that 


<1 (60) 


which relation again is true in practical cases. 


Eq. (58) now is changed into the following simpler 
form: 


P+ 


From Eq. (61), it follows that, as long as the cubic 
term is negligible, the relation between load and de- 
flection is of the form: 


E,=0 (61) 


(1 


This result is in agreement, as it should be, with explicit 
solutions of specific problems of the linear theory which 
have been obtained previously in reference 2. 


wom 
C2 E,h*t 


) 


+ cu(h/D? (t/P)(E/E + |= 0 (68) 


T?) (Ey/E.)(h ‘l)? (w 1)? 


~0O (59) 


1 + cs (ht/I?)(E,/G.) + (Ey/ 


The condition that the cubic term in Eq. (61) does 
not modify Eq. (62) to an appreciable extent is the fol- 
lowing: 


w < ce h/V1 + (63) 


Thus, when G, is of the same order of magnitude as 
E,, the load deflection relation is linear, as long as the 
deflection w is small compared with the thickness h of 
the composite plate. This is the same situation as for 
the homogeneous, isotropic plate. As the transverse 
shear deformability of the core layer becomes more 
pronounced, the range of linearity for w becomes 
smaller and smaller. We may ask what happens when 
G. — 0. We should expect that then the two face 
layers act as independent plates, and consequently the 
range of linearity becomes of order tf. However, we 
cannot carry out this transition G, — 0 in Eq. (63) 
because, in order to do this, we should have included in 
our analysis the transverse bending stiffness of the 
face layers themselves. We could, if we wished, ex- 
tend the present analysis in this direction so as to carry 
out the limiting process. For the present we content 
ourselves with the result, Eq. (63), which describes the 
effect of transverse core flexibility on the extent of the 
linear range of the load deflection characteristic, as long 


as the transverse bending rigidity of the individual 
face layer is negligible. 

As a further conclusion we have that in all but ex- 
treme cases we may omit the terms with e in the basic 
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on finite deflection theory as expressed in eqs. (40) to 
(44). 


5) THE Two SIMULTANEOUS EQUATIONS FOR AIRY’S 
SrrRESS FUNCTION F AND TRANSVERSE DEFLECTION Ww 


As long as we neglect the effect of transverse normal 
sress deformability in the present theory—that is, 
ete = O and E, © we may, following von Kar- 


| man's procedure for the homogeneous plate, arrive at 


asystem of two simultaneous equations as follows: 
We satisfy Eqs. (31) and (33) by setting 
N, = O°F/Oy’, N, = 0°F/0x? 
S = —0°?F/dxdy (64) 
The first of the two simultaneous equations is obtained 
by means of Eqs. (37) to (39), in the same form as for 
the homogeneous plate, 


= — (65) 


where V = 0?/Ox? + 07/Oy?. 

In order to obtain the second equation we first intro- 
duce Eqs. (32) and (34) into Eq. (35) and have, as for 
the homogeneous plate, 


YM, My 

t + N, = 0 (66) 


By means of the stress-strain relations, Eqs. (40) to 
(42), Eq. (66) is changed into 


, (Oa 
ag TP + 
Ox 


where D is the bending stiffness factor for the sandwich 
plate, 


D = (1/2)th*E,/(1 — vf) (68) 
According to Eqs. (43) and (44), we have 
—+—= — 69 
Ox hG, \ ox oy (69) 


We introduce Eq. (69) into Eq. (67) and observe once 
more Eq. (35). This gives, as the second of the two 
simultaneous equations, the following: 


dual 


ex- 
asic 


thE, 
DV?V = — ———, A?} X 
7 ¢ — ) 
Oy? Ox? Oxdy dxdy Ox? Oy? 


It is seen that the transverse shear deformability of the 
core introduces a group of new terms on the right side 
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differential eqs. (35) and (37) to (42), and there re- 
32)\f mains then the effect of transverse shear deformability 


of Eq. (70). When G, o Eq. (70) reverts to the 
known form of this equation for the homogenous plate. 

The solution of the system consisting of Eqs. (65) 
and (70) may be carried out by means of the procedures 
that have been developed for the case G, = © and will 


not be considered here. 


(6) OVERALL BUCKLING OF SANDWICH PLATES 


From Eq. (70) we may readily deduce an equation 
that governs the buckling of the composite plate. Let 
No, So, Ny be the components of the planar stress 
field, the stability of which is to be determined. From 
Eq. (70) the following differential equation for the 


_ buckling is derived: 
thE, 
ied ( 2(1 — ) 
O*w 
+235 N. 
( ‘ag (71) 


As an example, consider a simply supported rec- 
tangular plate subject to edge thrusts: 


Nu —No, Nin —M (72) 


Taking a plate of length a and width 0 and a deflection 
pattern 


(73) 
the following equation for the determination of the 
buckling loads is obtained : 


mr \* |? thE; 
l 
2 2 2 2 
a b a b 
In evaluating Eq. (74) we keep 6 fixed, assume JV, to 
be a specified positive or negative fraction of No, and 


determine that value of a which makes Np smallest. 
By treating the problem in this way we may also set 


w = sin (mmx/a) sin (may/b) 


m=n=1. We set as abbreviations 
th E; b? 
= 4 X = 75 
Eq. (74) can then be written in the form, 
Nob? (1 + A)? 
= = 76 


We further consider the two special cases, Ni/No = 
1 and Ni/No = 0. For Ni/No 1 we find that Vp is 
smallest when \ = 0 and that 


= 1/(1 + (77) 
For a square plate {\ = 1) we find 
= 2/(1 + (78) 


For Ni/No = 0 we find that, as long as uw < 1, No is 
smallest when 


. 
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A= (1+ — 


which means that transverse shear deformations. re- 
duce the wave length for the smallest buckling load. It 
is 


fla + w)/(1 u)] 4/(1 + »)*, 1 
(79) 


When yu 2 1, we have f = fin, for \ = © with the 
value 


= f() 1/u, 1 (79a) 


Ic is interesting to note that for this case of loading ff 
present theory becomes inadequate as soon as u becom 
of the order of magnitude unity, since it then giveg 
wave length for the buckling of the composite plq 
which eventually becomes of the order of the face lay 
thickness. 
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